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Abstract 

We introduce a new class of integrators for stiff ODEs as well as SDEs. Examples 
of subclasses of systems that we treat are ODEs and SDEs that are sums of two 
terms, one of which has large coefficients. These integrators are (i) Multiscale: they 
are based on flow averaging and so do not fully resolve the fast variables and have a 
computational cost determined by slow variables (ii) Versatile: the method is based 
on averaging the flows of the given dynamical system (which may have hidden slow 
and fast processes) instead of averaging the instantaneous drift of assumed sepa- 
rated slow and fast processes. This bypasses the need for identifying explicitly (or 
numerically) the slow or fast variables (iii) Nonintrusive: A pre-existing numerical 
scheme resolving the microscopic time scale can be used as a black box and easily 
turned into one of the integrators in this paper by turning the large coefficients on 
over a microscopic timescale and off during a mesoscopic timescale (iv) Convergent 
over two scales: strongly over slow processes and in the sense of measures over fast 
ones. We introduce the related notion of two-scale flow convergence and analyze 
the convergence of these integrators under the induced topology (v) Structure pre- 
serving: They inherit the structure preserving properties of the legacy integrators 
from which they are derived. Therefore, for stiff Hamiltonian systems (possibly 
on manifolds), they can be made to be symplectic, time-reversible, and symmetry 
preserving (symmetries are group actions that leave the system invariant) in all 
variables. They are explicit and applicable to arbitrary stiff potentials (that need 
not be quadratic). Their application to the Fermi- Pasta-Ulam problems shows ac- 
curacy and stability over four orders of magnitude of time scales. For stiff Langevin 
equations, they are symmetry preserving, time-reversible and Boltzmann-Gibbs re- 
versible, quasi-symplectic on all variables and conformally symplectic with isotropic 
friction. 
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1 Overview of the integrator on ODEs 

Consider the following ODE on M d , 



u 



G(u e 



+ -F(u e 
e 



(1.1) 



In Subsections 1.9, 2.1, 3.1 3.5 and 



4.1| we will consider more general ODEs, stiff de- 



terministic Hamiltonian systems (2.1), SDEs (( |3.l[ ) and (3.15)) and Langevin equations 
((4.1) and (4.2)); however for the sake of clarity, we will start the description of our 
method with (1.1). 

Condition 1.1. Assume that there exists a diffeomorphism r\ := (r) x ,rj y ), from ~R d onto 
M. d ~ p x W (with uniformly bounded C l ,C 2 derivatives), separating slow and fast vari- 
ables, i.e., such that (for all e > 0) the process (rcf,yf) = (rf (uf) , if (uf)) satisfies an 
ODE system of the form 

\x e = g(x e ,y e ) x e = x 
\y e = ]f(x e ,y e ) yo = yo 



(1.2) 



Condition 1.2. Assume that the fast variables in (1.2) are locally ergodic with respect 



to a family of measures /i drifted by slow variables. More precisely, we assume that there 
exists a family of probability measures fi(x, dy) on W indexed by x £ W i ~ p and a positive 
function T \— > E(T) such that lim^oo E(T) = and such that for all xo,yo,T and (j> 
uniformly bounded and Lipschitz, the solution to 



Y t = f(x ,Y t ) Y = y 



(1.3) 



satisfies 



1 



4>{Y s )ds- cP(y) f ,(x ,dy) < x(\\(xo,yo)\\)E(T)(\\<f>\\ L ^ + ||V0|U=o) (1.4) 
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where r i— >■ x(r) * s bounded on compact sets. 



Under conditions 1 1 . 1 1 and 1.2 it is known (we refer for instance to |95j or to Theorem 
14, Section 3 of Chapter II of |104j or to [88 j ) that x e converges towards xt defined as 
the solution to the ODE 



g(x,y)fj,(x,dy), x\ t=0 = x 



where fi(x, dy) is the ergodic measure associated with the solution to the ODE 

y = f(x,y) 



(1.5) 



(1.6) 



It follows that the slow behavior of solutions of (1.1) can be simulated over coarse time 



steps by first identifying the slow process x e and then using numerical approximations 



of solutions of (1.2) to approximate x e . Two classes of integrators have been founded on 
this observation: The equation free method |64| [65] and the Heterogeneous Multiscale 
Method [36, 40"tl35[l5]. One shared characteristic of the original form of those integrators 



is, after identification of the slow variables, to use a micro-solver to approximate the 



effective drift in (1.5) by averaging the instantaneous drift g with respect to numerical 



solutions of (1.6) over a time span larger than the mixing time of the solution to (1.6). 



1.1 FLAVORS 

In this paper, we propose a new method based on the averaging of the instantaneous flow 



of the ODE (1.1) with hidden slow and fast variables instead of the instantaneous drift 



of x e in ODE (1.2) with separated slow and fast variables. We have called the resulting 



class of numerical integrators FLow Averaging integratORS (FLAVORS). Since 



FLAVORS are directly applied to (1.1), hidden slow variables do not need to be identi- 



fied, either explicitly or numerically. F urth ermore FLAVORS can be implemented using 

i 

an arbitrary legacy integrator for (1.1) in which the parameter - can be controlled 
(figure [l]) . More precisely, assume that there exists a constant ho > such that $^ 
satisfies for all h < ho min(-, 1) and u G M. d 



u 



hG(u) - ahF(u)\ < Ch 2 (l + a) 



then FLAVOR can be defined as the algorithm simulating the process 
ut = ($£_ T o$)) fc (u ) for kS<t<(k + l)6 



(1.7) 



(1. 



where r is a fine time step resolving the fast time scale (r <C e) and 5 is a mesoscopic 
time step independent of the fast time scale satisfying r <C e <C 5 -C 1 and 



(-) 2 < 5 « 



(1.9) 



In our numerical experiments, we have used the "rule of thumb" 5 ~ where 7 is a 
small parameter (0.1 for instance). 
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Figure 1: A pre-existing numerical scheme resolving the microscopic time scale can be used as a black 
box and turned into a FLAVOR by simply turning on and off stiff parameters over a microscopic timescale 
t (on) and a mesoscopic timescale 5 (off). The bottom line of the approach is to (repeatedly) compose 
an accurate, short-time integration of the complete set of equations with an accurate, intermediate-time 
integration of the non-stiff part of the system. While the integration over short time intervals is accurate 
(in a strong sense), this is extended to intermediate time integration (in the sense of measures) using 
the interplay between the short time integration and the mesoscopic integration. The computational 
cost remains bounded independently from the stiff parameter 1/e because: (i) The whole system is only 
integrated over extremely short (r <C e) time intervals during every intermediate (<5) time intervals, (ii) 
The intermediate time step 5 (that of the non-stiff part of the system) is limited not by the fast time 
scales (e) but by the slow ones (0(1)). 



By switching stiff parameters FLAVOR approximates the flow of (1.1) over a coarse 
time step h (resolving the slow time scale) by the flow 



;i.io) 



where M is a positive integer corresponding to the number of "sa mple s" used to average 
the flow (5 has to be identified with -fe). We refer to subsection 



1.4 



for the distinction 



between macro and meso-steps, for the rationale and mechanism behind FLAVORS and 



the limits (1.9) 



Since FLAVORS are obtained by flow-composition, we will show in Section [2] and 
[4] that they inherit the structure preserving properties (for instance symplecticity and 
symmetries under a group action) of the legacy integrator for Hamlitonian systems and 
Langevin equations. 

is strongly accurate with 



Under conditions (1.9) on r and 5, we show that (1 



respect to (hidden) slow variables and weakly (in the sense of measures) accurate with 
respect to (hidden) fast variables . Motivated by this observation, we introduce the 
related notion of two-scale flow convergence in analogy with homogenization theory 
for elliptic PDEs [86, [3] and call it F-convergence for short. F-convergence is close in 
spirit to the Young measure approach to computing slowly advancing fast oscillations 
introduced in [TOl [9"]. 
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1.2 Two-scale flow convergence 

Let (Ct)teK+ be a sequence of processes on ~R d (functions from M + to M. d ) indexed by 
e > 0. Let (Xi)j gK + be a process on M d ~ p (p > 0). Let x \-> u(x, dz) be a function from 
M rf_p into the space of probability measures on M. d . 

Definition 1.1. We say that the process £f F-converges to v{Xt,dz) as e 1 and 
write £| )• v(Xt, dz) if and only if for all functions (p bounded and uniformly Lipshitz- 

e-S-0 

continuous on M d , and for all i > 0, 

r-t+h 



limlimi/ <p(£)ds= [ <p(z)v(X t ,dz) (1.11) 



1.3 Asymptotic convergence result 

Our convergence theorem requires that u| and do not blow up as e 1 0; more precisely, 
we will assume that the following conditions are satisfied. 

Condition 1.3. Assume that: 

1. F and G are Lipschitz continuous. 

2. For all uq, T > ; the trajectories (uf)o<t<r uniformly bounded in e. 

3. For all uq, T > 0, ifte trajectories (u~t)o<t<T «7"e uniformly bounded in e, < 5 < 
ho, t < min(roe, 8). 

For 7T, an arbitrary measure on M. d , we define n^ 1 * it to be the push forward of the 
measure 7r by rj^ 1 . 



Theorem 1.1. Let be the solution to (1.1) and ut be defined by (1.8). Assume that 
equation (1.1) and conditions \l . 1\ \lj\ and \i~^ are satisfied, then 



• u\ F-converges to n 1 * {Sx t ® M(-^*> dy)) as e 1 where Xt is the solution to 

X t = f g(X t ,y) fi(X t , dy) X = x . (1.12) 

• ut F-converges to rj^ 1 * {8x t ® fJ,(Xt,dy)) for e < 5/(— Cln<5), - i 0, -5 10 and 

(i) 2 Ho. 

Remark 1.1. The F-convergence of to ry" 1 * (<5x t <S> n(X t , dy)) can be restated as 

limlimi [ t+H tp(u e s )ds= f tp{rr\X t ,v))n{X t ,dy) (1.13) 
for all functions 99 bounded and uniformly Lipshitz-continuous on ~K d , and for all t > 0. 
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Remark 1.2. Observe that g comes from (1.5). It is not explicitly known and does not 



need to be explicitly known for the implementation of the proposed method. 

Remark 1.3. The limits on e, r and 5 are in essence stating that FLAVOR is accurate 



provided that r <C e (r resolves the stiffness of (1.1)) and equation (1.9) is satisfied. 

Remark 1.4. Throughout this paper, C will refer to an appropriately large enough con- 
stant independent from e, 5, r. To simplify the presentation of our results, we use the 
same letter C for expressions such as 2Ce c instead of writing it as a new constant C\ 
independent from e, 5, r. 



1.4 Rationale and mechanism behind FLAVORS 

We will now explain the rationale and mechanism behind FLAVORS. We refer to Sub- 
section 7.1 of the appendix for the detailed proof of Theorem |1.1| Let us sta rt b y 



considering the case where rj is the identity diffeomorphism. Let ip^ be the flow of (1.2). 
Observe that ip° (obtained from (p^ by setting the parameter - to zero) is the flow of 



(1.2) with y e frozen, i.e., 



<p°(x,y) = (xt,y) where x t solves ^ = g(x,y), x = x. (1-14) 



The main effect of FLAVORS is to average the flow of (1.2) with respect to fast degrees 

of freedom via splitting and re-synchronization. By splitting, we refer to the substitution 

i ' i 

of the flow ipt by composition of p®_ T and (p£ , and by re-synchronization we refer to 
the distinct time-steps 5 and r whose effects are to advance the internal clock of fast 

variables by r every step of length 5. By averaging, we refer to the fact that FLAVORS 

i 

approximates the flow (p^ by the flow 

p h :={<p\_ r o4) M (1-15) 

where h is a coarse time step resolving the slow time scale associated with x £ , M is a 
positive integer corresponding to the number of samples used to average the flow (5 is 
identified with -h) and r is a fine time step resolving the fast time scale, of the order of 
e, and associated with y € . In general, analytical formulae are not available for f° and 
<pe and numerical approximations are used instead. 

Observe that when FLAVORS are applied to systems with explicitly separated slow 
and fast processes, they lead to integrators that are locally in the neighborhood of 
those obtained with HMM (or equation free) methods with a reinitialization of the fast 
variables at macrotime n by their final value at macrotime step n — 1 and with only one 
microstep per macrostep [37], [39] • 

We will now consider the situation where r/ is not the identity diffeomorphism and 
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give the rationale behind the limits (1.9). 



$4 



'/ 



V 



O — T 



nS+r 



As illustrated in the above diagram, since (xt,yt 

ating the d 

(x n 5,y n s) ■ 



>,y)(n+l)<5 

rj(ut), simulating u n s defined in 



(1.8) is equivalent to simulating the discrete process 

*L°*r)>0,yo) 



(1.16) 



where 



n-=v°n°v- 1 (i-i7) 

Observe that the accuracy (in the topology induced by F-convergence) of u% with respect 



to ul, solution of ( 1.1 ), is equivalent to that of (xt, yt) with respect to (xf , yf) defined by 

(1.18) 



(1.2). Now, for the clarity of the presentation, assume that 

= u + hG{u) + a/iF(u) 



Using Taylor's theorem and (1.18), we obtain that 



*l{x,y) = (x,y) + h(g(x,y),0)+ah(0,f(x,y))+ [ v T Ressr,{u+tv)v(l -tf dt (1.19) 

J o 



with 



u := rj x (x, y) and 



h(G + aF)o7] 1 (x,y) 



1.20) 



It follows from equations (1.19) and (1.20) that ^ is a first order accurate integrator 



approximating the flow of (|1.2l) and 'i'lTis a first order accurate integrator approximating 



the flow of (1.14). Let h be a coarse time step and 5 a mesostep. Since x remains nearly 



constant over the coarse time step, the switching (on and off) of the stiff parameter ^ 
averages the drift g of x with respect to the trajectory of y over h. Since the coarse 
step h is composed of j mesosteps, the internal clock of the fast process is advanced by 
k x ~. Since h is of the order of one, the trajectory of y is mixing with respect to the 
local ergodic measure fj, provided that j- 3> 1, i.e. 



5<C 



;i.2i) 



Equation ( 1.21[ ) corresponds to the right hand side of equation (1.9). If rj is a non- linear 

diffeomorphism (with non-zero Hessian), it also follows from equations (1.19) and (1.20) 

^ 

that each invocation of the integrator ^ occasions an error (on the accuracy of the 

i 

slow process) proportional to (-) 2 . Since during the coarse time step h, ^ is solicited 
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I -times, it follows that the error accumulation during h is | x (^-) 2 . Hence, the accuracy 



of the integrator requires that I x (£) <C 1, i.e. 



' T \ 2 



(1.22) 



Equation (1.22) corresponds to the left hand side of equation (1.9). 



Observe that if rj is linear, its Hessian is null and the remainder in the right hand 



side of (1.19) is zero. It follows that if rj is linear, the error accumulation due to fine 



time steps on slow variables is zero and condition (1.21) is sufficient for the accuracy of 
the integrator. 

It has been observed in [38] and in Section 5 of [111] that slow variables do not need 
to be identified with HMM/averaging type integrators if the relation between original 
and slow variables is linear or a permutation and if 

At 



T 

M e 



(1.23) 



where is M the number of fine-step iterations used by HMM to compute the average the 
drift of slow variables and At is the coarse time step (in HMM) along the direction of 



the averaged drift. The analysis of FLAVORS associated with equation (1.19) reaches a 



similar conclusion if rj is linear in the sense that the error caused by the Hessian of r\ in 
(1.19) is zero and in the (sufficient) condition (1.21) is analogous to (1.23) for M = 1. It 
is also stated on Page 2 of [38] that "there are counterexamples showing that algorithms 
of the same spirit do not work for deterministic ODEs with separated time scales if the 
slow variables are not explicitly identified and made use of. But in the present context, 
the slow variables are linear functions of the original variables, and this is the reason why 
the seamless algorithm works. " Here, the analysis of FLAVORS associated with equation 



(1.19) shows an algorithm based on an averaging principle would indeed, in general, not 



work if rj is nonlinear (and (1.22) not satisfied) due to the error accumulation (on slow 
variables) associated with the Hessian of n. However, the above analysis also shows 



that if condition (1.22) is satisfied, then, although rj may be nonlinear, flow averaging 



integrators will always work without identifying slow variables. 



1.5 Non asymptotic convergence result 

Theorem 1.2. Under assumptions and notations of theorem there exists C > 
such that for 5 < ho, t < tq€ and t > 0, 

\x e t -rf{u t )\ < Ce ct Xl {uQ,e,5,T) 



and 



Tj t 



t+T 



tp(u s )ds- I ip(rj 1 (X t ,y))n{X t ,dy) 
Jrp 

< X2(uQ,e,5,T,T,t){\\ip\\ L °° + ||V^[|i«.) 
where xi and X2 are functions converging towards zero as e < <5/(Cln I] 



(1.24) 



ant 



d^YliO (andT 10 for X2 J- 



(1.25) 



9 



Remark 1.5. For e < 6/(— ClnS) and 5— + - < 1, the following holds 



(1.26) 



and X2 satisfies a similar inequality. 

Remark 1.6. Choosing r ~ je and (5 ~ 7-, where 7 is a small constant independent from 



e, Theorem 1.2 shows that the approximation error of FLAVOR is bounded by a function 
of 7 converging towards zero as 7 1 0. If follows that the speed up is of the order of 
r ~ ?, i.e., scales like ^ at fixed accuracy. In order to be able to compare FLAVOR with 
integrators resolving all the fine time steps we have limited the speed up in the numerical 
experiments to 200 x (but the latter can be arbitrary large as e ! 0). For sufficiently 
small e, we observe that FLAVORS with microstep r and mesostep 5 overperform their 
associated legacy integrator with the same microstep r over large simulation times (we 



refer to Section 6.3 on the Fermi-Pasta- Ulam problem). This phenomenon is caused 
by an error accumulation at each tick (microstep) of the clock of fast variables. Since 
FLAVORS (indirectly, i.e., without identifying fast variables) slow down the speed of 
this clock from 1 to a value ~ - independent from e this error does not blow up 
as e i (as opposed to an integrator resolving all fine time steps). Because of this 
reason, if this error accumulation on fast variables is exponential, then the speed up at 
fixed accuracy does not scale like -, but like where T is the total simulation time. 
A consequence of this phenomenon can be seen in Figure [HM (associated with the FPU 
problem) where Velocity Verlet fails to capture the 0(e~ 1 ) dynamics with a time step 
h = 10~ 5 whereas FLAVORS remain accurate with r = 10 -4 and 5 = 2- 10 -3 . 

Remark 1.7. The reader should not be surprised by the presence of the exponential 



factor e ct in (1.24). It is known that global errors for numerical approximations of 
ODEs grow, in general, exponentially with time (see for instance [55]). These bounds 
are, however, already tight; consider, for instance, how error propagates in a generic 
numerical scheme applied to the special system of x = x. It is possible to show that 
the increase of global errors is linear in time only for a restricted class of ODEs (using 
techniques from Lyapunov's theory of stability |115] ). Notice that the constant C in the 
exponential of our bound does not scale with e _1 , and therefore the bound is uniform 
and rather tight. 

Remark 1.8. We refer to [40J for higher order averaging based methods. In particular, 
|40j shows how, after identification of slow variables, balancing the different error con- 
tributions yields an explicit stable integration method having the order of the macro 
scheme. 



1.6 Natural FLAVORS 



Although convenient, it is not necessary to use legacy integrators to obtain FLAVORS. 
More precisely, theorems |1.1| and 1.2 remain valid if FLAVORS are defined to be algo- 
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rithms simulating the discrete process 



ut := (6f_ T o6 e T ) k {u ) for k5 < t < (k + 1)5 (1.27) 

where 6^. and 0f_ T are two mappings from M. d onto M. d (the former approximating the flow 
of the whole system (1.1) for time r, and the latter approximating the flow of v = G{v) 
for time 5 — r), satisfying the following conditions. 

Condition 1.4. Assume that: 

1. There exists ho,C > such that for h < ho and any u G M. d , 

\6%(u)-u-hG(u)\<Ch 2 (1.28) 



2. There exists tq, C > 0, such that for - < tq and any u G M d , 



6%{u) -u-tG{u) - -F(u) <C(-Y 



(1.29) 



3. For alluQ, T > 0, the discrete trajectories yyOs- 
bounded in e, < 5 < /iq, t < min(roe, 5). 



0<k<T/5 



are uniformly 



Observe that (1.8) is a particular case of (1.27) in which 9 e 

so 



9 G is obtained from the legacy integrator $ c 



and the mapping 



ay setting a to zero. 



1.7 Related work 

Dynamical systems with multiple time scales pose a major problem in simulations because 
the small time steps required for stable integration of the fast motions lead to large 
numbers of time steps required for the observation of slow degrees of freedom [109| 154] . 
Traditionally, stiff dynamical systems have been separated into two classes with distinct 
integrators: stiff systems with fast transients and stiff systems with rapid oscillations 
[HI E^l IHE] . The former has been solved using implicit schemes EH EH EH] , Chebyshev 
methods [701 [1] or the projective integrator approach [38]. These latter have been solved 
using filtering techniques [66J 1100] or Poincare map techniques [17J [90]. We also 
refer to methods based on highly oscillatory quadrature [32], [60J EH] , an area that has 
undergone significant developments in the last few years ]61j . It has been observed 
that at the present time, there exists no unified strategy for dealing with both classes 
of problems |35j . When slow variables can be identified, effective equations can be 
obtained by averaging the instantaneous drift driving those slow variables |104j . Two 
classes of numerical methods have been built on this observation: The equation-free 
method [BH [65] and the Heterogeneous Multiscale Method [361 HOI E3 E] ■ Observe that 
FLAVORS apply in a unified way to both stiff systems with fast transients and stiff 
systems with rapid oscillations, with or without noise, with a mesoscopic integration 
time step chosen independently from the stiffness. 
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1.8 Limitations of the method 



The proof of the accuracy of the method (theorems |1 , 1| and 1.2 ) is based on an averaging 



principle; hence, if e is not small (the stiffness of the ODE is weak), although the method 
may be stable, there is no guarantee of accuracy. More precisely, the global error of the 
method is an increasing function of e, 5, — , (^) 2 <5. Writing 7 := - the accuracy 

method requires 7 2 <C 5 -C 7. Choosing 5 = 72, the condition e <C 5 >C 1 (related to 

2 

computational gain) requires es <C 7 1 which can be satisfied only if e is small. 

The other limitation of the method lies in the fact that a stiff parameter - needs to 
be clearly identified. In many examples of interest (Navier-Stokes equations, Maxwell's 
equations,...), stiffness is a result of nonlinearity, initial conditions or boundary condi- 
tions and not of the existence of a large parameter 1 . Molecular dynamics can also create 
widely separated time-scales from non-linear effects; we refer, for instance, to [116J and 
references therein. 

1.9 Generic stiff ODEs 

FLAVORS have a natural generalization to systems of the form 

u^ = F(u a ' e ,a,e) (1.30) 
where u 1— > F(u, a, e) is Lipshitz continuous. 
Condition 1.5. Assume that: 

1. 6 1 — y F(u,a,e) is uniformly continuous in the neighborhood o/O. 

2. There exists a diffeomorphism rj := (n x ,r] y ) ! from M. d onto x MP, indepen- 
dent from e,a, with uniformly bounded C 1 , C 2 derivatives, such that the process 
(a#,i/£) = (?f (u"'°),ff satisfies, for all a>l, the ODE 

x Q = g(x a ,y a ) x% = x , (1.31) 

where g(x,y) is Lipschitz continuous in x and y on bounded sets. 

3. There exists a family of probability measures //(cc, dy) on MP such that for all 
Xo,yo,T ((xo,yo) := j?(uq)) and ip uniformly bounded and Lipschitz 







T M)ds- [ <p(y)n{xo,dy) <x(ll^o,yo)||)(^i(T) + ^ 2 (Ta 1 '))||V^|| L ^ 

JM.P 

(1.32) 

where r h-» x( r ) i s bounded on compact sets and E2 (r) — )• as r — ^ 00 and E\{r) — >• 
as r — > 0. 

4. For all uq, T > 0, the trajectories (u"'°)o<t<T are uniformly bounded in a > 1. 
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Remark 1.9. Observe that slow variables are not kept frozen in equation (1.32). The 



error on local invariant measures induced by the (slow) drift of x a is controlled by E2. 
More precisely, the convergence of the right hand side of ( |1 .32 ) towards zero requires the 
convergence of T towards zero and (at the same time) the divergence of Ta v towards 
infinity. 

Assume that we are given a mapping $^' e from ~R d onto M. d approximating the flow 
of (1.30). If the parameter a can be controlled then <3?^' e can be used as a black box for 



accelerating the computation of solutions of (1.30). 

Condition 1.6. Assume that: 

1. There exists a constant ho > such that <3? a ' e satisfies for all h < /iomin(^j, 1) 
< e < 1 < a 



|*h' e («) -u- hF(u, a,e)\< C{u)h 2 (l + a 2 ") 
where C(u) is bounded on compact sets. 



(1.33) 



2. For all Uq, T > 0, the discrete trajectories ^(^^ T ^"^(ito)^ are uni- 

formly bounded in < e < 1, < <5 < ho, r < min(/toe I/ , 5). 
FLAVOR can be defined as the algorithm given by the process 



Ut 



^ T o^'Y(uo) for k5 < t < (k + 1)5 



,e\k 



(1.34) 



The theorem below shows the accuracy of FLAVORS for 5 <C ho, r <C e u and (^7) <C 



Theorem 1.3. Let u^' e be the solution to (1.30) with a = 1/e and ut be defined by 



1.34). Assume that Conditions\l.5\and\l.6\are satisfied then 



• uf e F -converges towards 77 1 * (dx t ® fi(X t ,dy)) as e I where X t is the solution 
to 

X t = f g(X t ,y)ii(Xt,dy) X = x , (1.35) 

J MP 

• As e i 0, re - " i 0, i 0, -Jfj i 0, u t F -converges towards r/ _1 * (<5x t ®n{Xt, dy)) 
as e 1 where Xt is the solution of (1.35). 



Proof. The proof of Theorem 1.3 is similar to that of Theorem 1.1 and 3.1 Only the idea 
of the proof will be given here. The condition e <C 1 is needed for the approximation 
of u a,e by u a, ° and for the F-convergence of u^'°. Since yf = 7] y (u"'°) the condition 



r <C is used along with equation (1.33) for the accuracy of $r in (locally) approxi- 
mating yf. The condition 5 <C allows for the averaging of 5 to take place prior to a 

significant change of xf ; more precisely, it allows for m 3> 1 iterations of <&t ' prior to 
a significant change of . The condition (^7) <C 5 is required in order to control the 

error accumulated by m iterations of <3?t € - Q 
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2 Deterministic mechanical systems: Hamiltonian equa- 
tions 



Since averaging with FLAVORS is obtained by flow composition, FLAVORS have an 
inherent extension to multiscale structure preserving integrators for stiff Hamiltonian 
systems, i.e. ODEs of the form 

p=-d q H(p,q) q = d p H(p,q) (2.1) 

where the Hamiltonian 

H(q,p):=±p T M- 1 p+V(q) + ~U(q) (2.2) 

represents the total energy of a mechanical system with Euclidean phase space M. d x M. d 
or a cotangent bundle T*A4 of a configuration manifold A4. 

Structure preserving numerical methods for Hamiltonian systems have been devel- 
oped in the framework of geometric numerical integration [54\ 172] and variational inte- 
grators |78|, [75] . The subject of geometric numerical integration deals with numerical 
integrators that preserve geometric properties of the flow of a differential equation, and it 
explains how structure preservation leads to an improved long-time behavior [53]. Vari- 
ational integration theory derives integrators for mechanical systems from discrete vari- 
ational principles and are characterized by a discrete Noether theorem. These methods 
have excellent energy behavior over long integration runs because they are symplectic, 
i.e., by backward error analysis, they simulate a nearby mechanical system instead of 
nearby differential equations. Furthermore, statistical properties of the dynamics such 
as Poincare sections are well preserved even with large time steps |16j . Preservation of 
structures is especially important for long time simulations. Consider integrations of a 
harmonic oscillator, for example: no matter how small a time step is used, the ampli- 
tude given by Forward Euler / Backward Euler will increase / decrease unboundedly, 
whereas the amplitude given by Variational Euler (also known as symplectic Euler) will 
be oscillatory with a variance controlled by the step length. 

These long term behaviors of structure-preserving numerical integrators motivated 
their extension to multiscale or stiff Hamiltonian systems. We refer to [31] for a recent 
review on numerical integrators for highly oscillatory Hamiltonian systems. Symplectic 
integrators are natural for the integration of Hamiltonian systems since they reproduce at 
the discrete level an important geometric property of the exact flow |20] . For symplectic 
integrators primarily for (but not limited to) stiff quadratic potentials, we refer to the 
Impulse Method, the Mollified Impulse Method, and their variations |51|, 11101 196] . 
which require an explicit form of the flow map of stiff process. In the context of vari- 
ational integrators, by defining a discrete Lagrangian with an explicit trapezoidal ap- 
proximation of the soft potential and a midpoint approximation for the fast potential, 
a symplectic (IMEX — IMplicit-EXplicit) scheme for stiff Hamiltonian systems has been 
proposed in [105J. The resulting scheme is explicit for quadratic potentials and implicit 
for non quadratic stiff potentials. We also refer to Le Bris and Legoll's (Hamilton-Jacobi 
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derived) homogenization method |20j. Asynchronous Variational Integrators [73] pro- 
vide a way to derive conservative symplectic integrators for PDEs where the solution 
advances non-uniformly in time; however, stiff potentials require a fine time step dis- 
cretization over the whole time evolution. In addition, multiple time-step methods [106J 
evaluate forces to different extends of accuracies by approximating less important forces 
via Taylor expansions, but it has issues on long time behavior, stability and accuracy, as 
described in Section 5 of [73 J . Fixman froze the fastest bond oscillations in polymers to 
remove stiffness by adding a log term resemblant of entropy-based free energy to com- 
pensate [33]. This approach is successful in studying statistics of the system, but does 
not always reconstruct the correct dynamics |91| [59l [TI] . 

Several approaches to the homogenization of Hamiltonian systems (in analogy with 
classical homogenization [Til E2]) have been proposed. We refer to .M-convergence 
introduced in |101|, 115]. to the two-scale expansion of solutions of the Hamilton- Jacobi 
form of Newton's equations with stiff quadratic potentials |20j and to PDE methods in 
weak KAM theory (3JJ. We also refer to [25], [35] and [93] . 

Obtaining explicit symplectic integrators for Hamiltonian systems with non-quadratic 
stiff potentials is known to be an important and nontrivial problem. By using Verlet/leap- 
frog macro-solvers, methods that are symplectic on slow variables (when those variables 
can be identified) have been proposed in the framework of HMM (the Heterogeneous 
Multiscale Method) in (1021 124] . A "reversible averaging" method has been proposed 
in [TTJ for mechanical systems with separated fast and slow variables. More recently, 
a reversible multiscale integration method for mechanical systems was proposed in [6J 
in the context of HMM. By tracking slow variables, [6J enforces reversibility in all vari- 
ables as an optimization constraint at each coarse step when minimizing the distance 
between the effective drift obtained from the micro-solver (in the context of HMM) and 
the drift of the macro-solver. We are also refer to [99J for HMM symmetric methods for 
mechanical systems with a stiff potentials of the form - Y^j=\ 9j(l) 2 - 

2.1 FLAVORS for mechanical systems on manifolds 



Assume that we are given a first order accurate legacy integrator for (2.1) in which the 
parameter 1/e can be controlled, i.e. a mapping <I>^ acting on the phase space such that 
for h < ho min(l, a - 2 ) 

n(q,P)-(q,P)-HM- 1 p,-V(q)-aU(q)) <Ch 2 (l + a) (2.3) 



Write 0,5, the FLAVOR discrete mapping approximating solutions of (2.1) over time 
steps (5>e, i.e. 

(Q(n+l)5,P(n+l)s) == ©<5(<7n,5, Pud) ■ (2-4) 

FLAVOR can then be defined by 

e s := $g_ T o $| (2.5) 



Theorem 1.3 establishes the accuracy of this integrator under Conditions 1.5 and |1.6 



provided that r <C yfe <C 5 and ^- <C 5 -C 
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2.1.1 Structure preserving properties of FLAVORS 

We will now show that FLAVORS inherit the structure preserving properties of their 
legacy integrators. 

Theorem 2.1. If for all h, e > <& e h is symmetric under a group action, then Q$ is 
symmetric under the same group action. 

Theorem 2.2. // <I>^ is symplectic on the co-tangent bundle T*Ai of a configuration 



manifold Ai, then 0$ defined by (2.5) is symplectic on the co-tangent bundle T*A4. 

Theorem |2.1| and Theorem |2 ,2| can be resolved by noting that "the overall method is 
symplectic - as a composition of symplectic transformations, and it is symmetric - as a 
symmetric composition of symmetric steps" (see Chapter XIII. 1.3 of |54J). 

Write 

n := (<&_ fc ) -1 (2.6) 
Let us recall the following definition corresponding to definition 1.4 of the Chapter V of 



Definition 2.1. A numerical one-step method ^ is called time-reversible if it satisfies 
$1 = $h. 

The following theorem, whose proof is straightforward, shows how to derive a "sym- 
plectic and symmetric and time-reversible" FLAVOR from a symplectic legacy integrator 
and its adjoint. Since this derivation applies to manifolds, it also leads to structure- 
preserving FLAVORS for constrained mechanical systems. 

Theorem 2.3. If <J>^ is symplectic on the co-tangent bundle T*A4 of a configuration 
manifold M., then 

Q s ■- o $°'* T o $°_ T o $J (2.7) 

2 ~a- ~2~ 2 

is symplectic and time-reversible on the co-tangent bundle T*A4. 

Remark 2.1. Observe that (except for the first and last steps) iterating @$ defined by 



(2.7) is equivalent to iterating 



@S ■■= ®° S '* T ° $°s-t o $| o $|'* (2. 

2 2 



It follows that a symplectic, symmetric and reversible FLAVOR can be obtained in a 



nonintrusive way from a Stormer/Verlet integrator for (2.1) |53 |, I55 | |114| . 



2.1.2 An example of a symplectic FLAVOR 



If the phase space is M rf x then an example of symplectic FLAVOR is obtained from 



Theorem 2.2 by choosing to be the symplectic Euler (also known as Variational Euler 



or VE for short) integrator defined by 

XPJ \ -V(q)-aU(q) 



and letting ©5 be defined by (2.5) 
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2.1.3 An example of a symplectic and time-reversible FLAVOR 



If the phase space is the Euclidean space IR^xIR^, then an example of symplectic and time- 
reversible FLAVOR is obtained by letting O^ be defined by Equation (2.7) of Theorem 



2.3 



by choosing <J>^ to be the symplectic Euler integrator defined by (2.9) and 



+ h 



M' l p 

-V(q + hM-ip) - all(q + hM^p) 



(2.10) 



2.1.4 An artificial FLAVOR 



There is not a unique way of averaging the flows of (2.2 ). We present below an alternative 



method based on the freezing and unfreezing of degrees of freedom associated with fast 
potentials. We have called this method "artificial" because it is intrusive. With this 



method, the discrete flow approximating solutions of (2.1) is given by (2.4) with 



tr 



(2.11) 



where 9j is a symplectic map corresponding to the flow of H slow (q,p) := V(q), approx- 
imating the effects of the soft potential on momentum over the mesoscopic time step 5 
and defined by 

0Y{q,p) = [q,p-SVV(q)). (2.12) 

9% is a symplectic map approximating the flow of H^ ast (q,p) := hp 1 M~ 1 p + -U(q) over 
a microscopic time step r: 

9 € T {q,p) = {q + TM~ l p : p- T -WU{q + tM~ l p)) (2.13) 

&gL T is a map approximating the flow of the Hamiltonian H* ree (q,p) := ^p T M~ 1 p under 
holonomic constraints imposing the freezing of stiff variables. Velocities along the direc- 
tion of constraints have to be stored and set to be before the constrained dynamics, 
i.e., frozen, and the stored velocities should be restored after the constrained dynamics, 
i.e., unfrozen; geometrically speaking, one projects to the constrained sub-symplectic 
manifold, runs the constrained dynamics, and lifts back to the original full space. Of- 
tentimes, the exact solution to the constrained dynamics can be found (examples given 



in Subsections 5.3, 5.2 6.2 6.3 and 6.4) 



When the exact solution to the constrained dynamics cannot be easily found, one may 
want to employ integrators for constrained dynamics such as SHAKE [94J or RATTLE 
[3] instead. This has to be done with caution, because symplecticity of the translational 
flow may be lost. The composition of projection onto the constrained manifold (freez- 
ing), evolution on the constrained manifold, and lifting from it to the unconstrained 
space (unfreezing) preserves symplecticity in the unconstrained space only if the evolu- 
tion on the constrained manifold preserves the inherited symplectic form. A numerical 
integration preserves the discrete symplectic form on the constrained manifold, but not 
necessarily the projected continuous symplectic form. 
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Remark 2.2. This artificial FLAVOR is locally a perturbation of nonintrusive FLAVORS. 
By splitting theory [8T|l54"]. 



itr 



(2.14) 



whereas o $r 



- M.' r ^s^t ° ^ > where 0^ ree is the flow of H^ ree (q,p) under no 

constraint. The only difference is that constraints are treated in 6 tr but not in 9^. 

i 

Remark 2.3. This artificial FLAVOR can be formally regarded as . In contrast 



natural FLAVOR is $° 5 _ 



The advantage of this artificial FLAVOR lies in the fact that only r <C yfi <C 5 

2 

and 5 <C ^ are required for its accuracy (and not j < 5). We also observe that, 
in general, artifical FLAVOR overperforms nonintrusive FLAVOR in FPU long time 



(0(cj )) simulations (we refer to Subsection 6.3) 



2.2 Variational derivation of FLAVORS 

FLAVORS based on variational legacy integrators [78] are variational too. Recall that 
discrete Lagrangian is an approximation of the integral of the continuous Lagrangian 
over one time step, and Discrete Euler-Lagrangian equation (DEL) corresponds to the 
critical point of the discrete action, which is a sum of the approximated integrals. The 
following diagram commutes: 



Singlescale L 



FLAVORization 



action principle 



Multiscale 

action principle 



Singlescale DEL FLAVORlzaUon > Multiscale DEL 



For example, recall Variational Euler (i.e. symplectic Euler) for system (2.2) with 
time step h 

ipk+1 =Pk-h[VV(q k ) + lVU(q k )} 
1 q k +i = q k + hp k+ i 



(2.15) 



can be obtained by applying variational principle to the following discrete Lagrangian 

2 



Ld)!\qk,qk+i) = h 



1 / gfc+i - gfc 

2 I h 



V(q k ) + -U(q K 



(2.16) 



Meanwhile, FLAVORized Variational Euler with smallstep r and mesostep 5 



= Pk-r[VV(q k ) + \VU(q k )] 
= Qk + rp' k 
Pk+1 =p' k -(5-r)VV(q' k ) 
Qk+i =q' k + (S- r)p k+1 



(2.17) 
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can be obtained by applying variational principle to the FLAVORized discrete La- 
grangian 



Lds(qk,q' k , Qk+i) = L d l /e (q k , q' k ) + L d ° s _ T (q k , q k+1 ) 



I (q' h -q k 



V(q k ) + -U(q k ) 



+ (5-r) 



1 / gfc+i - q'k 

2 V S-t 



V{q' k 



(2.18) 



FLAVORizations of other variational integrators such as Velocity Verlet follow sim- 
ilarly. 



3 SDEs 

Asymptotic problems for stochastic differential equations arose and were solved simulta- 
neously with the very beginnings of the theory of such equations [ 104] . Here, we refer 
to the early work of Gikhman [39] , Krylov [67J EH] , Bogolyubov [13] and Papanicolaou- 
Kohler [87]. We refer in particular to Skorokhod's detailed monograph [104] . As for 
ODEs, effective equations for stiff SDEs can be obtained by averaging the instantaneous 
coefficients (drift and the diffusivity matrix squared) with respect to the fast compo- 
nents; we refer to Chapter II, Section 3 of |104j for a detailed analysis including error 
bounds. Numerical methods such as HMM [37j and equation-free methods [7J have been 
extended to SDEs based on this averaging principle. Implicit methods in general fail to 
capture the effective dynamics of the slow time scale because they cannot correctly cap- 
ture non-Dirac invariant distributions [76] (we refer to non-Dirac invariant distribution 
as a measure of probability on the configuration space whose support is not limited to 
a single point). Another idea is to treat fast variables by conditioning; here, we refer 
to optimal prediction [28] [271 [29] that has also been used for model reduction. We also 
refer to [SI E21 EDS EH E31 EE [2] . 

Since FLAVORS are obtained via flow averaging, they have a natural extension to 
SDEs developed in this section. As for ODEs, FLAVORS are directly applied to SDEs 
with mixed (hidden) slow and fast variables without prior (analytical or numerical) 
identification of slow variables. Furthermore, they can be implemented using a pre- 
existing scheme by turning on and off the stiff parameters. 

For the sake of clarity, we will start the description of with the following SDE on M. d : 

du\ = {G(u\) + -F(u\)) dt + {H(u\) + 4=#(«t)) dW t , ul = u (3.1) 
e V e 

where (Wt)t>o is a ci-dimensional Brownian Motion; F and G are vector fields on M d ; H 
and K ared x d matrix fields on M d . In Subsection 3.5, we will consider the more general 



form (3.15) 



Condition 3.1. Assume that: 

1. F,G,H and K are uniformly bounded and Lipschitz continuous. 
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2. There exists a diffeomorphism r\ := (rj x ,rj y ), from ~R. d onto ~R. d ~ p x W, independent 
of e, with uniformly bounded C , C 2 and C 3 derivatives, such that the process 
(xl,yf) = {r] x [uf) , r\ y {uf}) satisfies the SDE 

= g(x e ,y e )dt + a(x e ,y e )dW t , x e = x 
--y(x*,y*)dt + ±Q(x*,y e )dWt, y e = yo 

where g is d — p dimensional vector field; f a p-dimensional vector field; a is a 
(d — p) x d-dimensional matrix field; Q a p x d-dimensional matrix field and Wt a 
d-dimensional Brownian Motion. 

3. Let Yt be the solution to 

dY t = f(x ,Y t )dt + Q(x ,Y t )dWt Y = y (3.3) 

there exists a family of probability measures /j,(x, dy) on W indexed by x 6 W i ~ p 
and a positive function T i— > E(T) such that limr_ >00 E(T) = and such that for 
all xo,yo,T and <f> with uniformly bounded C r derivatives for r < 3, 



1 



E[0(y s )] - / 4>{y)n{x ,dy) < X (|[(^0, yo)\\)E(T) max \\4>\\ C r (3.4) 



where r i— > x( r ) * s bounded on compact sets. 



4- For all uq, T > 0, sup 0<t<T E 



is uniformly bounded in e. 



Remark 3.1. As in the proof of Theorem 1.1 the uniform regularity of F, G, H and K 



can be relaxed to local regularity by adding a control on the rate of escape of the process 
towards infinity. To simplify the presentation, we will use the global uniform regularity. 

We will now extend the definition of two-scale flow convergence introduced in Sub- 



section 1.2 to stochastic processes. 



3.1 Two-scale flow convergence for SDEs 

Let (£ e (i, w)) tgR+ wgfi be a sequence of stochastic processes on W 1 (progressively measur- 
able mappings from M + x Q to M d ) indexed by e > 0. Let (X t ) teR + be a (progressively 
measurable) stochastic process on M d_p (p > 0). Let x i— > u(x,dz) be a function from 
jjcZ-P j^to the space of probability measures on Mr. 

Definition 3.1. We say that the process £| F-converges to v{Xt, dz) as e \. and write 

£5 > v(Xt,dz) if and only if for all function ip bounded and uniformly Lipshitz- 

continuous on M d , and for all t > 

lim lim \ f + E ds = E[[ <p(z)u(X t , dz)] (3.5) 

h->0 e^Oh J t J R d 
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3.2 Non intrusive FLAVORS for SDEs 



Let w be a random sample from a probability spa ce ( SI, J 7 , IP) and $^(,w) a random 
mapping from M. d onto M rf approximating the flow of (3.1 ) for a = 1/e. If the parameter a 
can be controlled, then <3?^ can be used as a black box for accelerating the computation of 
solutions of (3.1 ) without prior identification of slow variables. Indeed, assume that there 



exists a constant ho > and a normal random vector £(u;) such that for h < ho min(^, 1) 



E 



| §l(u, oj)-u-hG(u)-ahF(u) - VhH(u)£(u) - VahK(u)£(cu) I 



< CTtlfl+oOi 



then FLAVOR can be defined as the algorithm simulating the stochastic process 



u = 

i 

ut = for k5 < t < (k + 1)5 



(3.6) 



(3.7) 



where u^, uj' k are i.i.d. samples from the probability space (O, J 7 , P), J < ho and r 6 (0, 5) 
such that r < Toe. Theorem |3.1| establishes the asymptotic accuracy of FLAVOR for 

t <§C e 5 and 

(-) f «5«-. 

3.3 Convergence theorem 



(3.8) 



Theorem 3.1. Let u e be the solution to (3.1) and ut defined by (3.1). Assume that 



equation (3.6) and Conditions \3.l\ are satisfied, then 

• u\ F -converges towards n~ l * (5x t &> n(Xt, dy)) as e 1 where Xt is the solution to 



dX t 



g(Xt,y)n(Xt,dy)dt + a(Xt)dB t X = x 



where a is a (d — p) x (d — p) matrix field defined by 

aa T = I aa T (x,y) ^(x,dy) 



(3.9) 



(3.10) 



and Bt a (d — p)- dimensional Brownian Motion. 

u t F -converges towards r/^ 1 * (6x t ® /J,(Xt,dy)) as e I 0, r < 5, | I 0, ^ I and 

(i) f Ho. 



The proof of convergence of SDEs of type (3.2) is classical, and a comprehensive 



monograph can be found in Chapter II of [104] . A proof of (mean squared) convergence 



of HMM applied to (3.2) (separated slow and fast variables) with a = has been 
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obtained in |37j . A proof of (mean squared) convergence of the Equation- Free Method 
applied to (3.2) with <t/0 but independent of fast variables has been obtained in [50] . 



Theorem 3.1 proves the convergence in distribution of FLAVOR applied to SDE (3.1) 



with hidden slow and fast processes. One of the main difficulties of the proof of Theorem 



3.1 lies in the fact that we are not assuming that the noise on (hidden) slow variables 
is null or independent from fast variables. Without this assumption, x\ converges only 
weakly towards Xt, the convergence of u £ can only be weak and techniques for strong 



convergence can not be used. The proof of Theorem 3.1 relies on a powerful result by 



Skorokhod (Theorem 1 of Chapter II of |104j ) stating that the convergence in distribution 
of a sequence of stochastic processes is implied by the convergence of their generators. 



We refer to Subsection 7.2 of the appendix for the detailed proof of Theorem 3.1 



3.4 Natural FLAVORS 

As for ODEs, it is not necessary to use legacy integrators to obtain FLAVORS for SDEs. 



More precisely, Theorem 3.1 remains valid if FLAVORS are defined to be algorithms 



simulating the discrete process 



ut = Uks for k6 < t < (k + 1)5 



(3.11) 



where Uk,oj' k are i.i.d. samples from the probability space and d% _and 6f_ 

are two random mappings from M. d onto M. d satisfying following conditions 



3.2 



More 



precisely, 9*(.,lj) approximates in distribution the flow of (3.1) over time steps r <C e. 
0h(.,uj) approximates in distribution the flow of 



dv$ = G{v e t ) dt + H{vl) dW t 



(3.12) 



over time steps h <C 1. 

Condition 3.2. Assume that: 

1. There exists ho, C > and a d- dimensional centered Gaussian vector £(w) with 
identity covariance matrix such that for h < ho, 



E 



6h(u,Lo) -u- hG{u) - \/hH{u)£{u) 



< Ch2 



(3.13) 



2. There exists ro,C > and a d-dimensional centered Gaussian vector with 
identity covariance matrix such that for - < tq, 



E 



:(u,oj)-u- tG(u) - -F(u) - y/rH(u)Z(u) 



T e K(u)^f 



<C7(I)i 
(3.14) 
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3. For all u , T > 0, sup < n < T/(5 E x{\\u n S 



r < min(roe, 6), where u is defined by (3.11). 



is uniformly bounded in e, < 5 < ho, 



3.5 FLAVORS for generic stiff SDEs 

FLAVORS for stochastic systems have a natural generalization to SDEs on M. d of the 
form 

du a ' e = F{u a ' e , a, e) dt + K(u a ' e , a, e) dW t (3.15) 

where (Wt)t>o is a (i-dimensional Brownian Motion, F and K are Lipshitz continuous 
in u. 



Condition 3.3. Assume that: 

1. 7 h-> F(u, a, 7) arad 7 a, 7) are uniformly continuous in the neighborhood 
o/O. 

2. There exists a diffeomorphism r\ := (rj x ,rj y ), from ~R. d onto ~R. d ~ p x W, independent 
from e, a, with uniformly bounded C , C 2 and C 3 derivatives, and such that the 
stochastic process (xf,yf) = (rj x (u"'°) , T] y (uf'° )) satisfies for all a > 1 the SDE 



dx c 



g(x a ,y a )dt + a(x a ,y a )dW t x£ = x 



(3.16) 



where g is d — p dimensional vector field, a is a [d — p) x d-dimensional matrix 
field, g and a are uniformly bounded and Lipschitz continuous in x and y. 

3. There exists a family of probability measures fi(x,dy) on MP such that for all 
xo,yo,T ((xo,2/o) := f/(^o)) andip with uniformly bounded C r derivatives for r < 3, 



^J\[p{y«)] ds-J <p(y)n(x ,dy) 



< X{\\(x ,y )\\) {E 1 {T)+E 2 {Ta v )) max \\ip\\ C r 

"(3.17) 



where r 1— > x( r ) * s bounded on compact sets and E 2 {r) — > as r — > 00 and E\{r) 
as r — > 0. 



^. For all uq, T > 0, sup 0<t<T E 



x(IK'°l 



is uniformly bounded in a > 1. 



Remark 3.2. As in the proof of Theorem |1.1| the uniform regularity of g and a can 
be relaxed to local regularity by adding a control on the rate of escape of the process 
towards infinity. To simplify the presentation, we have use the global uniform regularity. 

Let ubea random sample from a probability space (fi, F, P) an d <i>"' £ (., oS) a random 
mapping from M. d onto ' 



approximating in distribution the flow of (3.15 ) over time steps 

OS, 6 



t c. If the parameter a can be controlled, then ' can be used as a black box for 
accelerating the computation of solutions of (3.15). The acceleration is obtained without 
prior identification of the slow variables. 



Condition 3.4. Assume that: 
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1. There exists ho, C > and a d- dimensional centered Gaussian vector £(u>) with 
identity covariance matrix such that for h<ho,0<e<l<a and h < 
h mm(i 1) 



E 



\$% £ (u) -u-hF(u,a,e) -Vh£(uj)K(u,a,e)\ 2 ) <Ch§(l + a%) (3.18) 



2. For all uq, T > 0, sup < n <x/5 



E 



X(\\UnS\ 



t < min(hoe u ,5), where u is defined by (3.19). 



is uniformly bounded in e, < 5 < ho, 



FLAVORS Let 5 < h and r G (0,6) such that r < r e u . We define FLAVORS as 
the class of algorithms simulating the stochastic process 1 1->- u% defined by 



no = u 

ut = Uks for kS < t < (k + 1)5 



(3.19) 



where 0Jk,ui' k are i.i.d. samples from the probability space (fi, J 7 , P). 

Remark 3.3. simulates the randomness of the increment of the Brownian Motion 
between times <5/c and 5k + r. simulates the randomness of the increment of the 
Brownian Motion between times 5k + r and 8(k + l). The independence of uik and u' k is 
reflection of the independence of the increments of a Brownian Motion. 

The following theorem shows that the flow averaging integrator is accurate with 
respect to F-convergence for r <C t v <C 5 and 



-) 2 «<5« 



(3.20) 



Theorem 3.2. Let be the solution to (3.15) with a = 1/e and ut be defined by 



(3.19). Assume that Conditions\3.3\ and\3.4\ are satisfied then 



uf e F -converges towards rj 1 * (6x t &> fJ,(Xf,dy)) as e ], where Xt is the solution 
to 

dX t = I g{X t ,y)n{Xt,dy) + a{Xt)dB t X = x (3.21) 



where a is a (d — p) x (d — p) matrix field defined by 

aa T = j o-a T (x,y) n(x,dy) 
and Bt a (d — p)- dimensional Brownian Motion. 



(3.22) 



As e I 0, re ^4 0, 5^- \, 0, (^) 2 j | 0, F -converges towards r] 1 * (8x t 
n(Xt,dy)) as e \. where Xt is the solution to (3.21). 
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Proof. The proof of Theorem 3.2 is similar to the proof of Theorem 3.1 The condition 
e <C 1 is needed for the approximation of u a,e by u a, ° and for the i^-convergence of u^'°. 



Since yf 



r)»(u t 



a,0\ 



the condition r -C is used along with Equation (3.18) for the 



accuracy of <3?t' in (locally) approximating yf. The condition 5 <C ^ allows for the 
averaging of g and a to take place prior to a significant change of xat; more precisely, it 

allows for m ^> 1 iterations of &r ' e prior to a significant change of xat- The condition 

3 

(-rr) 2 < <5 is required in order to control the error accumulated by m iterations of 



4 Stochastic mechanical systems: Langevin equations 

Since the foundational work of Bismut [12] , the field of stochastic geometric mechanics 
has grown in response to the demand for tools to analyze the structure of continuous 
and discrete mechanical systems with uncertainty pm E3 E21 [301 E21 E3 E3 EH1 E3 H3 
El EH]. Like their deterministic counterparts, these integrators are structure preserving 
in terms of statistical invariants. 

In this section, FLAVORS are developed to be structure preserving integrators for 
stiff stochastic mechanical systems, i.e., stiff Langevin equations of the form 

---^V{q)dt-\VU{q)dt-cpdt+^/W^C2dW t 



and of the form 



M~ l p 

-VV{q) dt - \VU(q) dt -^pdt+ ^/2/3~ 1 ^dW t 



(4.2) 



where c is a positive symmetric d x d matrix. 



Remark 4.1. Provided that hidden fast variables remain locally ergodic, one can also 
consider Hamiltonians with a mixture of both slow and fast noise and friction. For the 



sake of clarity, we have restricted our presentation to (4.1) and (4.2) 



Equations (4.1) and (4.2) model a mechanical system with Hamiltonian 



H(q,p) := \p T M- x p + V(q) + -U{q). (4.3) 

The phase space is the Euclidean space M rf x W 1 or a cotangent bundle T*M of a 
configuration manifold A4. 

Remark 4.2. If c is not constant and A4 is not the usual ~K d x WL d Euclidean space, one 
should use the Stratonovich integral instead of the Ito integral. 
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4.1 FLAVORS for stochastic mechanical systems on manifolds 



As in Section [2j we assume that we are given a mapping <3?^ acting on the phase space 

(4.4) 



such that for h < ho min(l, a 2 ) 

K(<liP) ~ («>P) " h(M- l p, -V{q) - aU{q)) | < Ch\l + a) 
Next, consider the following Ornstein-Uhlenbeck equations: 

dp = -acpdt + \fa\j'lfi- x (hdWt 



The stochastic flow of (4.5) is defined by the following stochastic evolution map: 

, f2 



f*t r 



(4.5) 



(4.6) 



Let 5 < ho and r £ (0,(5) such that r < To/^/a. FLAVOR for (4.1) can then be defined 
by 

f(<Zo,Po) = (<Zo,Po) 



and FLAVOR for ((421) can be defined by 



(4.7) 



(Qo,Po) = (Qo,Po) 

1 1 

{Q(k+i)S,P(k+i)s) = $°_ T $r o ^l SkS+T (qkS,Pks) 



(4.8) 



3.4 



Theorem |3.2| establishes the accuracy of these integrators under Conditions 3.3 and 

r 



provided that r < ^/e < <5 and (^=) 2 < 5 < 



4.2 Structure Preserving properties of FLAVORS for stochastic me- 
chanical systems on manifolds 

1 

First, observe that if and VPf are symmetric under a group action for all e > 0, then 
the resulting FLAVOR, as a symmetric composition of symmetric steps, is symmetric 
under the same group action (see comment below Theorem 2.3). 

Similarly, the following theorem shows that FLAVORS inherits structure-preserving 
properties from those associated with (the component approximating the Hamilto- 
nian part of the flow) . 

Theorem 4.1. 



If <3?^ is symplectic, then the FLAVORS defined by (4 -7) and (4-8) are quasi- 
symplectic as defined in Conditions RL1 and RL2 of \8Jfl fa degenerates to a 
symplectic method if friction is set equal to zero and the Jacobian of the flow map 
is independent of (q,p) ). 
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• If in addition c is isotropic then FLAVOR defined by (J^.l) is conformally symplec- 
tic, i.e., it preserves the precise symplectic area change associated to the flow of 
inertial Langevin processes fEUjj. 

Proof. Those properties are a consequence of the fact that FLAVORS are splitting 
schemes. The quasi-symplecticity and symplectic conformallity of GLA has been ob- 
tained in a similar way in \17\ . □ 



4.2.1 Example of quasi-symplectic FLAVORS 

An example of quasi-symplectic FLAVOR can be obtained by choosing <I>^ to be the sym- 



plectic Euler integrator defined by (2.9). This integrator is also conformally symplectic 
if c is isotropic and friction is slow. 



4.2.2 Example of quasi-symplectic and time-reversible FLAVORS 



Defining <3?^ by (2.9) and by (2.10), an example of quasi-symplectic and time- 



reversible FLAVOR can be obtained by using the symmetric Strang splitting: 



(q{k+i)8iP(k+i)s) 



H+f,(fc+l)5 



4)0 9 —ft- o 9 



' kS,kS+ 



s_(q,p) (4.9) 



for (4.1 ) and 



(g(*+l)i,P(fc+l)*) = % +1) s- h{k+1)S ° *r* ° *°£r ° $ V *| ° n s , kS+ r(q,p) (4.10) 



for (4.2). This integrator is also conformally symplectic if c is isotropic and friction is 
slow. 



4.2.3 Example of Boltzmann-Gibbs reversible Metropolis-adjusted FLA- 
VORS 

Since the probability density of ^ti,t 2 can De explicitly computed, it follows that the 



probability densities of (4.9) and (4.10) can be explicitly computed, and these algo- 



rithms can be metropolized and made reversible with respect to the Gibbs distribution 
as it has been shown in [19] for the Geometric Langevin Algorithm introduced in |17j . 
This metropolization leads to stochastically stable (and ergodic if the noise applied 
on momentum is not degenerate) algorithms. We refer to [19J for details. Observe 
that if the proposed move is rejected, the momentum has to be flipped and the accep- 
tance probability involves a momentum flip. It is proven in [19] that GLA [T7] remains 
strongly accurate after a metropolization involving local momentum flips. Whether this 
preservation of accuracy over trajectories transfers in a weak sense (in distributions) to 
FLAVORS remains to be investigated. 
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Stability domain of non-intrusive FLAVOR 




(a) Nonintrusive FLAVOR 



3 
2.5 
2 

» 1.5 
1 



Stability domain of artificial FLAVOR 











Stable 


Unstable 







(b) Artificial FLAVOR 



Figure 2: Stability domain of non-intrusive and artificial FLAVOR applied to (5.1 1 as a function of 5 
and r/e. u = l/y^ = 1000. 



5 Numerical analysis of FLAVOR based on Variational Eu- 
ler 

5.1 Stability 

Consider the following linear Hamiltonian system 



tt { \ 1 9 1 2 -1-2 ^ f \2 

H{x,y,p x ,py) = -p x + -p y + -x +y(y-x) 
with cj 3> 1. Here 2^ is the slow variable and y — x is the fast variable. 



(5.1) 



It can be shown that, when applied to (5.1), Symplectic Euler (2.9) is stable if and 



only if h < V2/w. Write &s,t the non-intrusive FLAVOR (2.5) obtained by using Sym- 
plectic Euler (|2.9l) as a Legacy integrator. Write @$ T the artificial FLAVOR described 
in Subsection 12. 1.41 

Theorem 5.1. The non-intrusive FLAVOR &s,t with 1/\/t S> cu 3> 1 is stable if and 
only ifSe (0,2). 

The artificial FLAVOR 6g T with 1/r > u > 1 is stable if and only if 5 £ (0, 2y/2). 
Proof. The numerical scheme associated with 0,5 T can be written as 



(5.2) 



Vn+l 




Vn 


Xn+1 


= T 


%n 


(Py)n+l 




(Py)n 


(Px)n+1_ 




_(Px)n_ 
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The characteristic polynomial of T is 

A 4 + (-4 + 5 2 - 5 2 t 2 + 25t 3 - r 4 + 25tuj 2 - 5 2 t 2 oj 2 + 25t 3 uj 2 - rV)A 3 + (6 - 25 2 
+ 25 2 t 2 - 45t 3 + 2r 4 - 4cW + 5 3 tuj 2 + 25 W - 4<5r 3 w 2 - 5 3 r 3 u; 2 + 2rV 
+ 2<S 2 rV - 5t 5 uj 2 )X 2 



+ (-4 + <5 2 



5 2 t 2 + 2<5r 3 - r 4 + 2<5tw 2 - <5 2 r 2 w 2 + 25r 3 u; 2 - r 4 w 2 )A + 1 



(5.3) 



Since oj 3> 1, r <C 1/uj 2 , as long as 5 < 1 roots to the above polynomial are close to 
roots to the asymptotic polynomial 



A 4 + (5 2 - 4)A 3 + (6 - 25 2 )A 2 + (<5 2 - 4)A + 1 



(5.4) 



which can be shown to be 1 with multiplicity 2 and ^(2 — S 2 ± 5\/5 2 — 4). It is easy to 
see that all roots are complex numbers with moduli less or equal to one if and only if 
1*1 < 2. 

The numerical scheme associated with @g T can be written as in (5.2) with 

5- 
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0" 




"1 





r 


0" 







1 













1 





T 
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(5.5) 



The characteristic polynomial of T is 

2A 4 +(4w 2 r 2 +T(5+(5 2 -8)A 3 +(12-25 2 -25r-8rV+25 2 rV)A 2 +(4w 2 T 2 +Tj+(5 2 -8)A+2 

(5.6) 

Since w > 1, t < 1/u, as long as S < 1 roots to the above polynomial are close to roots 
to the asymptotic polynomial 



2A 4 + (5 2 - 8)A 3 + (12 - 25 2 )X 2 + (5 2 - 8)A + 1 



(5.7) 



which can be shown to be 1 with multiplicity 2 and |(4 — 5 2 ± 8\/5 2 — 8). All roots are 
complex numbers with moduli less or equal to one if and only if | S\ < 2^2 □ 

Figures 2(a)| and 2(b) illustrate the domain of stability of nonintrusive FLAVOR 
(based on symplectic Euler (2.5) and (2.9)) and artificial FLAVOR (2.11) applied to the 
flow of (5.1 ), i.e. values of 5 and r/e ensuring stable numerical integrations. We observe 
that artificial FLAVOR has a much larger stability domain than nonintrusive FLAVOR. 
Specifically, for nonintrusive FLAVOR and large values of S, t = o(y/e) is not enough 
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and one needs r = o(e) for a stable integration, whereas artificial FLAVOR only requires 
r = y/2e, a minimum requirement for a stable symplectic Euler integration of the fast 
dynamics. 

Notice that there is no resonance behavior in terms of stability; everything below the 
two curves is stable and everything outside is not stable (plots not shown). 



5.2 Error analysis 



The flow of (|5.1l) has been explicitly computed and compared with solutions obtained 



from nonintrusive FLAVOR based on symplectic Euler ((2.5) and (2.9)) and with arti- 



ficial FLAVOR (2.11). 



The total simulation time is T = 10, and absolute errors on the slow variable have 
been computed with respect to the Euclidean norm of the difference in positions between 
analytical and numerical solutions. Stability is investigated using the same technique 
used in Subsection 5.1 Figures 3(a) and |3(b)| illustrate errors as functions of mesostep 8 
and renormalized small step r/e. Observe that given 8 errors are minimized at specific 
values of r/e for both integrators, but the accuracy of nonintrusive FLAVOR is less 
sensitive to r/e. Figures [3(c) and |3(d)] plot the optimal value of r/e as a function of 8 
and the associated to error. Observe also that for nonintrusive FLAVOR the dependence 
of the optimal value of r/e on 8 is weak, whereas for artificial FLAVOR the optimal value 
of r/e roughly scales linearly with 8. Figure 3(e) and |3(f )] describe how error changes 



with smallstep r for mesostep 8 fixed. Figure 3(e) can be viewed in correspondence 



with the condition 8 « r/e required for accuracy. This requirement, however, is just 
a sufficient condition to obtain an error bound, as we can see in Figure [3(f)] There 
the weak dependence of the error on r/e for a fixed 8 shows that one does not have to 
choose the microstep with too much care or optimize the integrator with respect to its 
value, if artificial FLAVOR is used. As a matter of fact, all the numerical experiments 



illustrated in this paper (except for Figures 3(c) and 3(d) ) have been performed without 
any tuning of the value r/e. We have simply used the rule of thumb 8 ~ 7^ where 7 is 
a small parameter (0.1 for instance). 

Therefore, it appears that the benefits of artificial FLAVORS lie in their superior 
accuracy and stability. 

Notice that there is no resonant value of 8 or r. 



5.3 Numerical error analysis for nonlinear systems 

In this subsection, we will consider the nonlinear Hamiltonian system 

H(x,y,z,p x ,p y ,p z ) = -pi + -pi + -pi + x 4 + e _1 y(y - xf + e _1 y(z - y) 2 (5.8) 

Thus, the potential is U = ^(y - x) 2 + f(z - y) 2 and V = x 4 . Here acts as a 

slow degree of freedom and y — x and z — y act as fast degrees of freedom. 

Figure [4] illustrates t 1— )■ X M + ^M +Z W ( s \ ow variable, convergent strongly) and t 1— > 
(y(t) — x(t), z(t) — y{t)) (fast variables, convergent in measure) computed with symplec- 
tic Euler and with the induced symplectic FLAVOR ( |2.5| )). Define q := (x,y,z). To 
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(a) Error of nonintrusive FLAVOR as a (b) Error of artificial FLAVOR as a function 

function of 5 and r/s/e. Notice that not all of S and r/^/e 
pairs of step lengths lead to stable integra- 
tions. 
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(c) Optimal t / \fe and error of nonintrusive (d) Optimal r/^/e and error of artificial 
FLAVOR as functions of S FLAVOR as functions of S 

Error as a function of t for 3=0.01 



Error as a function of i for 8=0.01 




0.05r 

0.04- 
0.03- 
0.02- 
0.01 - 




(e) Error dependence on r/^/e for a given S: (f) Error dependence on r/^/e for a given S: 



nonintrusive FLAVOR 



artificial FLAVOR 



Figure 3: Error analysis of (|5.1|. Parameters are lj = y/e = 10 3 , a;(0) = 0.8 and y(0) = a;(0) + 1.1/u 
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Slow variable integrated by VE 



Fast variables integrated by VE 



Is 




10 20 30 40 50 
Time 

Slow variable integrated by FLAVORS 
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0.05 



-0.05 




0.05 
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Figure 4: Comparison between trajectories integrated by Variational Euler and FLAVOR (defined by 



(2.51 and ( 2.9 1) . FLAVOR uses mesostep 8 = 0.01 and microstep r = 0.0005, and Symplectic Euler uses 



time step r = 0.0005. Time axes in the right column are zoomed in (by different ratios) to illustrate the 
fact that fast variables are captured in the sense of measure. FLAVOR accelerated the computation by 
roughly 20x ( 5 = 20r). In this experiment e = 10~ 6 , oji = 1.1, u 2 = 0.97, x(0) = 0.8, y(0) = 0.811, 
z(0) = 0.721, p x (0) = 0, p y (0) = and p z (0) = 0. Simulation time T = 50. 



illustrate the .F-convergence property of FLAVOR, we fix H = 1, vary the mesostep 
5 = H/M by changing M and show the Euclidean norm error of the difference between 



M 



Y^i=o 1 Q(T ~ ih/M) computed with FLAVOR and computed with symplectic Euler 
in Figure 5(a)| Notice that without an averaging over time length h, the error will be no 
longer monotonically but oscillatorily decreasing as 5 changes (plots not shown) , because 
fast variables are captured only in the sense of measure. As shown in Figure 5(a) the 
error scales linearly with jj for M not too small, and therefore the global error is a 
linear function of the mesostep 5 and the method is first order convergent. Figure 5(b)| 
shows that the error in general grows linearly with the total simulation time, and this 
linear growt h of the error has been observed for a simulation time larger than u (e -1 / 2 ). 
Figure 5(c) shows that the error does not depend on uj ( e" 1 / 2 ) for a fixed 5, as long as e 
is not too large (i.e. uj not too small); the error is instead controlled by M. This is not 
caused by reaching the limit of machine accuracy, it is a characteristic of the method: 
the plateau for large uj corresponds to the complete scale separation regime of FLAVOR 
as a multiscale method. 

Notice that there is no resonant value of 5 in the sense of convergence. 
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simulation time T 



Figure 5: Error dependence on parameters in a FLAVOR simulation of (5 



The fact that the error scales linearly with total simulation time is a much stronger 
(numerical) result than our (theoretical) error analysis for FLAVORS (in which the 
error is bounded by a term growing exponentially with the total simulation time). We 
conjecture that the linear growth of the error is a consequence of the fact that FLAVOR 
is symplectic and is only true for a subclass of systems, possibly integrable systems. A 
rigorous analysis of the effects of the structure preservation of FLAVORS on long term 
behavior remains to be done. 

6 Numerical experiments 

6.1 Hidden Van der Pol oscillator (ODE) 

Consider the following system ODEs 



Mr cos 9 + r sin# — lr 3 cos 3 9) cos 9 — er cos 9 sin# 



-e cos 



(cos 9 + sin l 



^r 2 cos 3 9) sin ( 



(6.1) 



where eCl. Taking the transformation from polar coordinates to Cartesian coordinates 
by [x,y] = [r s'm8, r cos#] as the local diffeomorphism, we obtained the hidden system: 



-ey 

](x + y-ly 3 ) 



(6.2) 



Taking the second time derivative of y, the system can also be written as the 2 nrf -order 
ODE: 



y + y 



l 



(i 



y 2 )y- 



3.3) 



The latter is a classical Van der Pol oscillator [113J. Nonintrusive FLAVOR as defined by 



(1.34) can be directly applied to (6.1) (with hidden slow and fast processes) by turning 



on and off the stiff parameter K More precisely, defining $ e ' a (r, 6) by 

(r cos 9 + r sin 6 — lr 3 cos 3 6) cos 



I r 3, 

(cos 9 + sin 9 — ir 2 cos 3 $) sin ( 



r cos sin \ . , 
fh{ cosH (6 - 4) 
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(a) Transformed System with Separated Timescales in (x,y): Forward Euler: Timestep=5e-005 




500 1000 1500 2000 2500 3000 3500 4000 4500 5000 

(b) Transformed System with Separated Timescales in (x,y): FLAVORS: Timestep=0.01 




500 1000 1500 2000 2500 3000 3500 4000 4500 5000 

(c) System with Hidden Separation of Timescales in (r,6): FLAVORS: Timestep=0.01 



r(t)cos6(t) 




500 1000 1500 2000 2500 3000 3500 4000 4500 5000 

Time 



Figure 6: Over a timespan of 5/e (a) Direct Forward Euler simulation of (6.2 1 with time steps resolving 
the fast time scale (b) (nonintrusive ( 1.34 1) FLAVOR simulation of (6.2 1 (c) Polar to cartesian image 
of the (nonintrusive ( |1.34[ )) FLAVOR simulation of ( |6.1| l with hidden slow and fast variables. Forward 
Euler uses time step h = 0.05e = 0.00005. The two FLAVORS simulations use 5 = 0.01 and r = 0.00005. 
Parameters are \ = 1000, x(0) = 1, y(0) = 1 



FLAVOR is defined by (1.34) with u := (f,0), i.e. 



(f t ,9 t ) = (^ T o^ e ) k (r ,e ) for k5 < t < (k + 1)8. (6.5) 

We refer to Figure [6] for a comparison of integrations by Forward Euler, used as a 
benchmark, and FLAVORS. FLAVORS gives trajectories close to Forward Euler and 
correctly captures the O(^) period |113| of the relaxation oscillation. Moreover, a 200x 
acceleration is achieved using FLAVOR. 

6.2 Hamiltonian system with nonlinear stiff and soft potentials 



In this subsection, we will apply the Symplectic Euler FLAVOR defined by (2.5) and 



(2.9) to the mechanical system whose Hamiltonian is 

H{y,x,p y ,p x ) := + -pi + e~V + (x - y) 4 (6.6) 
Here, stiff potential e _1 ?7 = e _1 y 6 and soft potential V = (x — y) 4 are both nonlinear. 
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Figure 7: In this experiment, e = 10~ e , y(0) = 1.1, x(0) = 2.2, p y (0) = and ^(0) = 0. Simulation 



time T = 2. FLAVOR (defined by (2.5 I and (2.9 1) uses mesostep 8 = 10 and microstep r = 10 



Variational Euler uses small time step r 



10" 



and IMEX uses mesostep 8 



ur 



Since the fast 



potential is nonlinear, IMEX is an implicit method and nonlinear equations have to be solved at every 
step, and IMEX turns out to be slower than Variational Euler. FLAVOR is strongly accurate with 
respect to slow variables and accurate in the sense of measures with respect to fast variables. Comparing 
to Symplectic Euler, FLAVOR accelerated the computation by roughly lOOx. 



Figure [7] illustrates t i-> y(t) (dominated by a fast process), t \-t x(t) — y(t) (a slow 
process modulated by a fast process), and t h-» H(t) computed with: Symplectic Euler, 



the induced symplectic FLAVOR and (2.9)), and IMEX [105]. Notice that x-y is 



not a purely slow variable but contains some fast component, and therefore the FLAVOR 
integration of it contains a modulation of local oscillations, which could be interpreted 
as that fast component slowed down by FLAVOR. It's not easy to find a purely slow 



variable or a purely fast variable in the form of ( 1.2 ) for this example, but the integrated 



trajectory for such a slow variable will not contain these slowed-down local oscillations. 

I 'ii i2 ■•■ <Wi q 2m 

I stiff soft 

Figure 8: Fermi-Pasta-Ulam problem ~ ID chain of alternatively connected harmonic stiff and 
non-harmonic soft springs 
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Expansions of stiff springs 
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(a) By Variational Euler with small (b) By artificial FLAVOR pTlip with 
time step r' = 5 x 1CP J = 0.05/cj. 38 pe- mesostep S — 0.002 and microstep r — 
riods in Subplot2 with zoomed-in time 10~ 4 = 0.1/cj. 38 periods in Subplot2 
axis (^380 in total over the whole sim- with zoomed-in time axis (^380 in total 
ulation span). over the whole simulation span). 



Figure 9: Simulations of the FPU problem over T = 2uj. Subplot2 of both figures have zoomed-in time 
axes so that whether phase lag or any other distortion of trajectory exists could be closely investigated. 
In this experiment m = 3, u = 10 3 , x(0) = [0.4642,-0.4202,0.0344,0.1371,0.0626,0.0810] is randomly 
chosen and y(Q) = [0, 0, 0, 0, 0, 0] . 



6.3 Fermi-Pasta-Ulam problem 

In this subsection, we will consider the Fermi-Pasta-Ulam (FPU) Problem [32] illustrated 
by Figure [8] and associated with the Hamiltonian 

^ m 2 m m 

H{q,p) := - J2(P2i-l+P 2 2i) + X " *X-l) 2 + X^ 2 ^ 1 ~ ^) 4 " ( 6 - 7 ) 

i=l i=l i=0 

The FPU problem is a well known benchmark problem |79[ [S^j for multiscale integrators 
because it exhibits different behaviors over widely separated timescales. The stiff springs 
nearly behave like harmonic oscillators with period ~ Then, the centers of 

masses linked by stiff springs (i.e., the midpoints of stiff springs) change over a timescale 
0(1). The third timescale, 0(cu), is associated with the rate of energy exchange between 
stiff springs. The fourth timescale, 0(uj 2 ), corresponds to the synchronization of energy 
exchange between stiff springs. On the other hand, the total energy of the stiff springs 
behaves almost like a constant. This wide separation of timescales can be seen in Figure 



[9j [TOj and [L2j where four subplots address different scales: Subplot 1 shows the fast 
variables (q2i~ <?2i-i)/ V^; Subplot2 shows one of the slow variables (92+91)/ V^; Subplot3 
shows the energy transfer pattern among stiff springs, which is even slower; Subplot4 
shows the near-constant total energy of three stiff springs. All four subplots are time- 
series. A comprehensive survey on FPU problem, including discussions on timescales 
and numerical recipes, can be found in |54j . 

t igures [9(a)] and [9(b)] compare symplectic Euler (with time steps fine enough to 



resolve FPU over the involved long time scale) and with the artificial FLAVOR (2.11) 



On a timescale 0(u>) (uj 3> 1), FLAVOR captured slow variable's periodic behavior with 
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the correct period and phase, as well as the slower process of energy transfer. At the 
same time, FLAVOR accelerated the computation by roughly 40x (since 5 = 40r'). 

It is not worrisome that artificial FLAVOR produces stiff spring energy trajectories 
with rapid local oscillations, which exhibit in both thicker individual energy curves 
and total energy with larger variance. In fact, these local oscillations do not seem 
to affect the global transfer pattern nor its period and are caused by the numerical 
error asociated with microstep r. This can be inferred by using the artificial FLAVOR 
introduced in Subsection 2.1.4 with Q% corresponding to the exact flow of H^ ast (rather 
than its Variational Euler approximation: this specific artificial Euler resembles the 
Impulse Method, but the Impulse Method will yield unbounded trajectories if one runs 
even longer time simulations, whereas FLAVORS do not seem to have an error growing 
exponentially with the total simulation time) . As illustrated in Figure 1 1 , exact flow 
helps to obtain thin energy curves of stiff springs with no rapid local oscillations as 
well as a total energy with a variance smaller than that given by fine Variational Euler 
(Figure 9(a)), with the transfer pattern similar to Figure [9(b)| 

Now, we reach further to 0(uj 2 ) total integration time to investigate different inte- 
grators' performances in capturing the synchronized energy exchange pattern (Figure 



10j). 

There is a significant difference among stiff spring energy transfer patterns produced 
by Velocity Verlet, FLAVOR, IMEX and the Impulse Method. Here, there is no analytic 
solution or provably accurate method for comparison. FLAVOR is the only method 
that shows periodic behavior on the long time scale and convergence tests show that 
FLAVOR'S trajectories remain stable under small variations of step sizes. Notice that 
mathematically it can be shown that the dynamical system admits periodic orbits (for 
example, by Poincare-Bendixson theorem). Furthermore, it is physically intuitive that 
the three stiff springs should alternatively obtain their maximal and minimal energies, 
and these maximal energies should be of fixed values. In addition, if we change the 
slow potential to be quadratic the system is still very similar to non-harmonic FPU but 
now analytically solvable. There, the energy exchanging pattern (Figure 12) resembles 
FLAVORS' result of the non-harmonic system but not the other integrators'. Notice that 
if run on the modified quadratic FPU problem however, FLAVORS, Velocity- Verlet, 
IMEX and the Impulse Method all obtain perfect results (plots omitted). These are 
numerical evidences supporting FLAVOR on the 0(lo 2 ) timescale. 

It is worth discussing why Velocity- Verlet with a time step much smaller than the 
characteristic length of the fast scale (C(l/w))is still not satisfactory. Being a second 
order method, it has an error bound of 0(e T h 2 ). On the other hand, backward error 
analysis guarantees that the energy of the integrated trajectory oscillates around the 
true conserved energy, hence eliminating the possibility of exponential growth of the 
numerical solution. Nevertheless, at this moment there is no result known to the authors 
to link these two analytical results to guarantee long term accuracy on the stiff springs' 
energies. The energy exchange among stiff springs is in fact an delicate phenomenon, 
and a slight distortion in stiff spring lengths could easily disrupt its period or even its 
periodicity. 
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(a) By Velocity Verlet with small time step (b) By artificial FLAVOR (iTTty with 
h = 10 -5 . mesostep 8 = 0.002 and microstep r = 0.0005 = 

0.1/w. 




(c) By IMEX with mesostep 8 — 0.002. (d) By Impulse Method with mesostep 8 = 

0.002. 



Figure 10: Simulations of FPU problem over T = \uj 2 . Initial conditions are x(0) = [1,0,0, 0,0] 
and j/(0) = [0, 0, 0, 0, 0, 0] so that energy starts concentrated on the leftmost soft and stiff springs, 
propagates to the right, bounces back, and oscillates among springs. We chose a smaller uj = 200 
because with a larger u) it would take weeks to run Velocity Verlet on a laptop. 



These numerical observations seem to indicate that symplectic FLAVORS may have 
special long time properties. Specifically, although we could not quantify the error 
here because there is no benchmark to compare to when the total simulation time is 
0{u 2 ), the long term behavior seems to indicate an error growing much slower than 
exponentially (please refer to Remark 1.7 for a discussion on exponential error bounds 
and Figure [5(b) for another example of conjectured linear error growth). A rigorous 
investigation on FLAVORS' long time behavior remains to be done. 

Figure 13 summarizes FLAVOR'S performance on various timescales in a comparison 
to Velocity Verlet. 

Notice that there are many sophisticated methods designed for integrating the FPU 
problem (see |54j for a review), as well as general multiscale methods that can be applied 
to the FPU problem. HMM as one state-of-art method in the latter category, together 
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Figure 11: By artificial FLAVORS (Subsection 



1000 2000 3000 4000 



Ml 

00 9000 10000 

2.1.4 1 based on exact fast flow with mesostep S = 

0.002 and microstep r = 10" 4 . Less oscillatory Figure 12: Harmonic FPU, T = 50a;, exact solu- 
stiff spring energies. 38 periods in Subplot2 with tion 
zoomed-in time axis (~380 in total over the whole 
simulation span). 



with identification of slow variables [5] captured the energy transfer between stiff springs 
over a time span of the order of oj. Simulations shown here are over a time span of the 
order of C(w 2 ). 

6.3.1 On resonances 

Multiscale in time integrators are usually plagued by two kinds of resonances. 

The first kind, called Takens resonance |1U7| , is related to the fact there are no closed 
equations on slow variables [15]. FLAVORS avoid Takens resonance because, thanks to 
-F-convergence, the information on the local invariant measure of fast variables is not 
lost. Observe that the FPU problem has Takens resonance (the eigenfrequencies of 
the strong potential are identical). Nevertheless, FLAVORS still capture the solution 
trajectories given any large value of oj with mesostep <5> 1/u independent of u. 

The second kind [25J is related to instabilities created by interactions between pa- 
rameters e, r and S. For instance, if e _1 = lu 2 resonance might happen at co5 or wt 
equal to multiples of tt/2. The analysis provided in Section [5] shows that such unstable 
interaction does not occur, either in the sense of stability or in terms of peaks of error 
function. This can be intuitively understood by observing that FLAVORS never approx- 
imate cos(6lj), while on the other hand, it does approximate cos(ro;) whose resonance 
frequency r = 2it/uj is ruled off by the requirement of r <C e for nonintrusive FLAVOR 
and r < y/e for artificial FLAVOR. 

6.4 Nonlinear 2D primitive molecular dynamics 

Now consider a two-dimensional, two degrees of freedom example in which a point mass 
is linked through a spring to a massless fixed hinge at the origin. While the spring as well 



39 



Expansions of stiff springs 



Expansions of stiff springs 



(a) FLAVOR 



0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0. 

Midpoint position of the first stiff spring 




50 100 150 200 250 300 350 400 450 500 

Energies of stiff springs 



0(1) 

o( w ) 



1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 



0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.0B 0.09 

Midpoint position of the first stiff spring 




1000 2000 3000 4000 5000 6000 7000 



(b) Velocity Verlet 



9000 10000 



Figure 13: Quantities of interest in integrations of FPU over different timescales. FLAVOR (2.111 
captures the fastest timescale in the sense of measure, while Velocity Verlet cannot accurately capture 
the slowest (0(uj 2 )) timescale despite of the small time step it uses. Here FLAVOR is 200 times faster 



than Velocity Verlet. All parameters are the same as in Figure 10(a) and 10(b) e.g. u) = 200, S = 0.002, 
r = 0.0005 and h = 10 -5 . 



as the point mass are allowed to rotate around the hinge (the spring remains straight), 
the more the spring-mass tilts away from its equilibrium angle the more restorative force 
it will experience. This example is a simplified version of prevailing molecular dynamics 
models, in which bond lengths and angles between neighboring bonds are both spring- 
like; other potential energy terms are ignored. 

Denote by x and y the Euclidean coordinates of the mass, and p x , p y corresponding 
momentums. Also introduce polar coordinates {r,9), with x = rcos9 and y = rsin#. 
Then the Hamiltonian reads 



H = \vl + \v\ + \j{r ~ r ) 2 + (cos#) 2 



= k + \pI + \^\Vx^ - r ) 2 + -^— 2 (6.8) 
1 I y 2 x z + y z 

where tq is equilibrium bond length parameter and oj is large number denoting bond 
oscillation frequency. 

Remark 6.1. This seemingly trivial example is not easy to integrate. 

1. If the system is viewed in Euclidean coordinates (x, y,p x ,Py) it is completely non- 
linear with a nonpolynomial potential, and hence the Impulse Method or its vari- 
ations [5H II 101 UH [96] , or IMEX [105] , or the homogenization method introduced 
in [20] cannot be applied using a mesostep. 



2. If the Hamiltonian is rewritten in generalized coordinates (r,9,p r ,pQ), H = + 

2 

+ \oj 2 {t — ro) 2 + I cos(6*) 2 , a fast quadratic potential can be identified. 
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However, the mass matrix 







is not constant, but rapidly oscillating, and 



hence methods that work for quasi-quadratic fast potentials (i.e. "harmonic oscil- 
lator" with a slowly changing frequency) (|20j for example) cannot be applied. 




Figure 14: Simulation of (6.8 1. Symplectic Euler uses small time step r = 0.0002 and the induced 
symplectic FLAVOR ((2.51 and ( 2.9 1) uses mesostep 8 = 0.01 and microstep r = 0.0002. In this 
simulation u = 500, x(0) = 1.1, t/(0) = 0.8, p x (0) = 0, p„(0) = and simulation time T = 100. 



Figure |14| compares symplectic Euler with the induced symplectic FLAVOR ((2.5) 



and (2.9)) applied to (6.8) in Euclidean coordinates. 

FLAVOR reproduced the slow 6 trajectory while accelerating the simulation time by 
roughly 50x (since 5 = 50r). It can also be seen from both energy fluctuations and the 
trajectory of the fast variable that the fast process' amplitude is well captured although 
its period has been lengthened. 



6.5 Nonlinear 2D molecular clipper 

We now consider a united-atom representation of a three atom polymer with two bonds 
(e.g. propane or water molecule). This is a simplified version of several prevailing 
molecular dynamics force fields (for example, CHARMM [21] . AMBER [33|, or a simpler 
example of butane \92\ 193]). Using conservation of momentum, we fix the coordinate 
system in the 2D plane defined by the three atoms. Introduce both Cartesian coordinates 
(xi,yi,x 2 ,y2,X3,y 3 ), as well as generalized coordinates r\ = y 7 (x 2 - xi) 2 + (y 2 - yi) 2 
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Figure 15: One exemplary configuration of a propane molecule 



and r 2 = 1/ (X3 — x 2 ) 2 + (1J3 — y 2 ) 2 for bond lengths and 6 for the angle between the two 
bonds (Figure 15). The kinetic energy is 



K.E. 



t -2 , -2\ , 1 / * 2 1 * 2 \ 1 r -2 , -2\ 

-m 1 (x 1 +y l ) + -m 2 (x 2 + y 2 ) + -m 3 (x 3 + y 3 ) 



3.9) 



where mi,m2,and 1113 denote the masses of the atoms. 

The potential energy consists of a bond term and a bond angle term, both of which 
are of harmonic oscillator type: 



RE. = Vh 



bond 



+ v c 



angle 



Vbond = 2 K r[(ri ~ r f + (r 2 - r ) 2 ] 

Vangle = ^K g (cOs(6) - COs(6*o)) 2 



(6.10) 
(6.11) 

(6.12) 



Notice that the system is in fact fully nonlinear: if written in generalized coordinates, 
the kinetic energy will correspond to a nonlinear and position dependent mass matrix, 
whereas in Cartesian coordinates, both terms in the potential energy are non-polynomial 
functions in positions. 

In the case of propane, mi = 15/i, m<i = 14/z, 7713 = 15/i where \x = 1.67 • 10 kg, 
r = 1.531, K r = 83.7kcal/(molA 2 ), 9 = 109.5° and K e = 43.1kcal/mol [92]. 

The propane system is characterized by a separation of timescales to some extent: 
bond stretching and bond-angle bending are characterized by 10 14 and 10 13 Hz vibra- 
tional frequencies respectively |117j . For investigation on FLAVORS, we use unitless 
parameters and exaggerate the timescale separation by setting K r to be 8370 and Kg to 
be 4.31. We also let [i = 1 without loss of generality for arithmetic considerations. 

In this system, the bond potential is the fast potential and the bond-angle potential is 
the slow one. It is well known that using a large time step at the timescale corresponding 
to the bond-angle potential by freezing bond lengths produces biased results, and many 
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physics based methods have been proposed to remedy this difficulty (for example by the 
approach of Fixman [33]; also see a review in |117j ). On the other hand, few multiscale 
methods work for this fully nonlinear system. 
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Figure 16: Simulations of exaggerated propane molecule (Subsection 



ler uses h — 0.01 and the induced symplectic FLAVOR ((2.51 and (2.9 ) parameters are 8 = 



6.51. Symplectic Eu- 



0.1 and r = 0.01. Initial conditions are [xi, J/i, X2, V2, %3, Jte] = [0,0,1.533,0,2.6136,1.0826] and 
[m 1 x 1 ,miy 1 ,m 2 x 2 , m 2 y 2 , m 3 x 3 , m 3 y 3 ] = [-0.4326, -1.6656, 0.1253, 0.2877, -1.1465, 1.1909]. 



Figure |l6| compares symplectic Euler with the induced symplectic FLAVOR ((2.5) 



and (2.9)) applied in Euclidean coordinates. lOx acceleration is achieved. A simula- 



tion movie is also available at http://www.cds.caltech.edu/~mtao/Propane.avi and 
http : //www Tacm . caltech . edu/~owhadi/| 



6.6 Forced nonautonomous mechanical system: Kapitza's inverted pen- 
dulum 

As the famous Kapitza's inverted pendulum shows [63J (for recent references sec |6j 
for numerical integration and [§7] for generalization to the stochastic setting), the up 
position of a single pendulum can be stabilized if the pivot of the pendulum experiences 
external forcing in the form of vertical oscillation. Specifically, if the position of the 
pivot is given by y = sin(ut), the system is governed by 

19= {g + u} 2 sm(2irut)}sm0 (6.13) 
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where denotes the clockwise angle of the pendulum from the positive y direction, I is 
the length of the pendulum and g is the gravitational constant. In this case, the rapid 
vibration causes the pendulum to oscillate slowly around the positive y direction with a 
0(1) frequency. 



VE + d'Alambert 



VE + d'Alambert 





Figure 17: Simulations of the inverted pendulum. The integration by Variational Euler + 
d'Alembert principle uses time step h = 0.2/uj/V~1 « 0.000067, while FLAVOR (defined 
by ( |6.15| » uses 5 = 0.002 and r = 0.2/u/Vl. Also, g = 9.8, I = 9, 0(0) = 0.2, 0(0) = 
and oj = 1000 

A single scale integration of this system could be done by Variational Euler with 
discrete d'Alembert principle for external forces jTS] 



(6.14) 



fi = oj 2 sin(27ra;i/i) 
Pi+i =Pi + h[g + fi] sin 0j 
9i+i = 0i + hp i+1 /l 

where the time step length h has to be smaller than 0(1/uj). 
FLAVOR is given by 

q n5+T = q nS + rpns/l 

PnS+r = PnS + rg sm.(q n5+T ) + u? sin(27rwnr) 

Q(n+1)8 = QnS+r + (S - T)p nS+T /l 

P(n+i)S = PnS+r + ($- r)gsm(q {n+1)S ) 



(6.15) 
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Observe that the time dependent force is synchronized on the r time scale instead of 
the 6 time scale, specifically uj 2 sin(27ro;nr) instead of u 2 sm(2irum5) in (6.15) 



Numerical results are illustrated in Figure 17 (also available as a movie at http:// 



www. cds . caltech.edu/~mtao/InvertedPendulum.avi and http : //www. acm. caltech 



edu/~owhadi/). Notice in this example that 9, being the only degree of freedom, con- 
tains a combination of slow and fast dynamics. FLAVOR could only capture the fast 
dynamics in the sense of measures, and this is why dents appear as modulation on the 
slow oscillation of 9. On the other hand, although this forced system does not admit a 
conserved energy, the value of the Hamiltonian should oscillate periodically due to the 
periodic external driving force. While a non mechanics based method such as Forward 
Euler often produces an unbounded growth or a decrease in the energy, FLAVORS do 
not have this drawback. 

Remark 6.2. Consider the case of a rapid potential of the form Q 2 (qi)q 2 /e 2 (where q\ is 
the slow an q2 the fast variable). In the limit of a vanishing e, it is known that the term 
contributes to the effective Hamiltonian with a contribution V{q\) (the so called Fixman 
term). On may have the intuition that FLAVOR would only be consistent with a term 
of the form 7V(<Zi), where < 7 < 1 is only a fraction of one, because the rapid force 
is only accounted for over a time r < 5. This intuition is not correct because the effect 
of FLAVOR is not to account the rapid force for over a time r < 5 but to slow down 
the rapid force by a fraction r/e. This effect can also be seen in the algorithm (6.15) 
where the force term to 2 sm(2irumT) has been slowed down by a factor t/S (the Kapitza's 
inverted pendulum illustrates a similar phenomenon where rapid oscillations contribute a 
stabilizing term to the effective Hamiltonian, nevertheless FLAVORS remain accurate). 



6.7 Nonautonomous SDE system with hidden slow variables 

Consider the following artificial nonautonomous SDE system 



du 
dv 



3(u+v) 2 
4 



-\ (^pf + 5sin(2vrt)) dt 

1 / v—u \2 
"2 \ 2 



+ 5 sin(2vrt) ) dt + - 



1 ( (u+v^3 
u+v \ 3 



+ c- 
+ c- 



v — u 
2 



v—u 
2 



dt- J-dWt 

dt + J\ dW t 

(6.16) 



3(u+v)' 2 

where c is a positive constant and the two dWt terms refer to the same Brownian Motion. 
The system (6.16) can be converted via the local diffeomorphism 

J u = (x — c) 1 / 3 



y 

I v = (x — c) 1 / 3 + y 
into the following hidden system separating slow and fast variables 



dx 



1„,2 



y 2 dt + 5sm(27rt)dW t 



dy = Ux-y)dt+ J 2 dW t 



(6.17) 



(6.18) 



Nonintrusive FLAVOR (3.7) can be directly applied to (6.16) using a time step 5 » e 
without prior identification of the slow and fast variables, i.e., without prior identification 
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(a) FLAVOR: the system in (u,v) with slow and fast variables hidden 




(c) Euler-Maruyama: the transformed system in (js,y) with an explicit separation of slow and fast variables 
2 I 1 i 1 1 1 1 1 1 1 




Time 



Figure 18: (a) Integration of ( |6.16[ ) by nonintrusive FLAVOR (3.71 using mesostep step S = 0.01 (b) 
Integration of (6.161 by Euler-Maruyama using fine time step h = 10~ 4 (c) Integration of (6.181 by 
Euler-Maruyama using the same small step h — 10~ 4 . Expectations of the slow variable (whether or not 
hidden) are obtained by empirically averaging over an ensemble of 100 independent sample trajectories, 
e = 10 -4 , x(0) = 1 + e, y(0) — 1, T = 2 (the expectation of the real solution will blow up around T = 3). 
We have chosen c = 10 so that the transformation is a diffeomorphism. 



of the slow variable x or of the system (6.18). The expected values of solutions of (6.16) 



integrated by FLAVORS with mesostep 5 and Euler-Maruyama with a small time step 



r are presented in Figure 18 FLAVOR has accelerated the computation by lOOx. 



6.8 Langevin equations with slow noise and friction 

In this subsection, we consider the one dimensional, two degrees of freedom system 
modeled by the SDEs (now both springs are quartic rather than harmonic): 



dy = pydt 
dx = p x dt 

dp y = -e~ 1 y 3 dt - A(y - xfdt - cp y dt + adW} 
dp x = — 4(x — y) 3 dt — cp x dt + adW? 



(6.19) 



We compare several autocorrelation functions and time-dependent moments of this 



stochastic process integrated by quasi-symplectic FLAVOR ((4.7) and (2.9)) and Geo- 



metric Langevin Algorithm (GLA) [TTj. FLAVOR and GLA gave results in agreement 
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Figure 19: SDE ( |6.19| ): autocorrelation functions of E[y(i)y(0)] (dominantly fast) and 
of E[(x(i) — y(t))(x(0) — y(0))] (dominantly slow), empirically obtained by GLA and 
FLAVORS. 



(Figure [l9j [20(a]| [20(b)] ). Since GLA is weakly-convergent and Boltzmann-Gibbs pre- 



serving, this is numerical evidence that quasi-symplectic FLAVOR is also. 

Expectations are empirically calculated by averaging over an ensemble of 100 sample 
trajectories with T = 30, e = 10" 8 , r = 0.001, 5 = 0.01. y(0) = 2.1/u (with u := 1/y/e), 
x(0) = y(0) + 1.8, c = 0.1 and a = 0.5. GLA uses time step h = 0.001. Noise and 
friction are slow here in the sense that they are not of the order 0(u) or larger. 

As shown in the plots, in the regime dominated by deterministic dynamics (roughly 
from t = to t = 8) various moments calculated empirically by FLAVORS and GLA 
are in agreement, indicating that the same rate of convergence towards the Boltzmann- 
Gibbs distribution is obtained. And in that regime, autocorrelation functions of the 
slow variables agree, serving as numerical evidence that FLAVORS is weakly converging 
towards the SDE solution, whereas autocorrelation functions of the fast variables agree 
only in the sense of measures (after time averaging over a mesoscopic (i(l)) time span). 
The fluctuations between FLAVORS and GLA for large time are an effect of the finite 
number of samples (100) used to compute sample averages. 

Recall that if the noise is applied to slow variables, FLAVORS do not converge 
strongly but only in the sense of distributions. 
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Figure 20: SDE (6.19): Empirical moments obtained from simulations of ensembles of 



(6.19) with GLA and quasi-symplectic FLAVOR (Subsection 4.2.1) 



6.9 Langevin equations with fast noise and friction 

Consider a system with the same configuration as above. The difference is that the soft 
spring oscillates at a frequency nonlinearly dependent on the stiff spring's length, and 
the left mass experiences strong friction and noise while the right mass does not. The 
Hamiltonian is 



H(y,x,p y ,p x ) 
and the governing SDEs are: 



1 



1 



1 



p I y + - P i/2 + -y y /i + ey(x-y) 



(6.20) 



dy = 
dx - 

dp y 

dpx 



■ Pydt 
-- p x dt 

= -u A y 3 dt -(2 + y 
= -2(x - y)e y dt 



x)(y — x)e y dt — uj 2 cp y dt + uoodW 1 



(6.21) 



In this system, the deterministic dynamics and the effects of noise and friction both 
involve a O(l/co 2 ) timescale. We have implemented the fast noise and friction version 



of FLAVORS ( fl4.8| and \2.9\). 

In Figure [2l| we have plotted the first and second moments of the dominantly slow 
variable x(t) — y(t) as well as the first moment of the dominantly fast variable y{t) 
as functions of time. Moments of the dominantly slow variable integrated by quasi- 



symplectic FLAVOR (Subsection 4.2.1) and GLA [T7] concur, numerically suggesting 
weak convergence and preservation of Boltzmann-Gibbs. lOOx computational accelera- 
tion is achieved. 
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GLA Empirical Expectation of Slow Variable GLA Empirical Expectation of Fast Variable GLA Empirical Variance of Slow Variable 




FLAVORS Empirical Expectation of Slow FLAVORS Empirical Expectation of Fast FLAVORS Empirical Variance of Slow 
2 i ■ 1 0.015 i ■ 1 4 




Figure 21: E[x( t) - y(t)], E[j/(t)], and E[x(t) - y(t)} 2 obtained by GLA and quasi-symplectic FLAVOR 
(Subsection 4.2.1 1. Expectations are empirically calculated by averaging over an ensemble of 50 sample 



trajectories with T = 10, u = 100, r = 10~ 4 , S = 0.01. y(0) = x(0) = y(0) + 1.8, c = 0.1 and 

a = 1. GLA uses time step h = 10 -4 . 



7 Appendix 

7.1 Proof of theorems 11.11 and 11.21 

Define the process t h-» (x t , yt) by 

{xt,yt) ■= v(u t ). 



(7.1) 



It follows from the regularity of r\ that it is sufficient to prove the F-convergence of 
(xt,yt) towards Sx t ® fi(X t ,dy). Moreover, it is also sufhcient to prove the following 
inequalities (7.2), (7.3) in order to obtain inequalities (1.24) and (1.25) 

\x e t -x t \ < Ce ct i)i (uo,e,8,T) 

and 



(7.2) 



t+T 



ip(x s ,y s )ds - / ip(X t ,y)n(X t ,dy) < ^ 2 (u , e, 5, r, T, t)(\\ip\\ L °° + \\ V<^||l«>) 
Jrp 

(7.3) 

Now define by 



(7.4) 
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Define if)? by 

i>h( x ,y) : =v o &h or f 1 {x,y) ( 7 - 5 ) 

Proposition 7.1. The vector fields f and g associated with the system of Equations 



1.2) are Lipschitz continuous. We also have 



(x t ,y t ) = {ips_ T °ipr) ( x O:Vo) for k5<t<(k + l)5. 
Moreover, there exists C > such that for h < ho and - < tq we have 



Ts2 



\r T (x, y) - (x, y) - r(g(x, y), 0) - - (0, f(x, y)) \ < C(-) 



(7.6) 
(7.7) 

t 

y) - (x, y) - h(g(x, y),0) \ < Ch 2 . (7.8) 

Furthermore, given xo,yo, the trajectories of (x|,y|) and (xt,y~t) <ire uniformly bounded 
in e, 5 < ho, t < min(roe, 5). 



anc 



Proof. Since (x, y) = rj(u), we have 

x = (G + -F)Vrf o rf 1 (x, y) 
e 

y = (G+ -F)Vrj y o V ~ l (x,y). 
Hence, we deduce from Equation ( |1.2[ ) of Condition 1 1 . 1 1 that 

g(x,y) = GVrf o rj~ l (x,y) 
f(x, y) = FVr? o r,~ l (x, y). 



(7.9) 
(7.10) 

(7.11) 
(7.12) 



We deduce the regularity of / and g from the regularity of G, F and n. Equation (7.6) 
is a direct consequence of the definition of and i^ 9 h and Equation (1.27) (we write 
(^0)2/0 ) := ^(^o))- Observe that Equation (1.2) of Condition |1 . 1| also requires that 

FVrf = GVrf = 0. (7.13) 

Now observe that 

y) - ( x , v) - {a(x, y),o)r- (o, f(x, y)) - = 

(7 1 oe t T - v -T(GVr ! x ,0) - ^(0,FVr] y )) oif\x,y). 
Using ( 7.13| ), (1.29), Taylor expansion and the regularity of rj we obtain (7.7). Similarly 
^(x, y) - (x, y) - h[g(x, y),0) := (tjo0°- V (x, y) - h[GVrf, 0)) o rf l (x, y). (7.15) 



(7.14) 



Using (7.13), (1.28), Taylor expansion and the regularity of r\ we obtain ( |7.8| ). The 
uniform bound (depending on xo^yo) on the trajectories of {x\,yf) and (xtiVt) is a 
consequence of the uniform bound (given uq) on the trajectories of u\ and Uf D 
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It follows from Proposition 7.1 that it is sufficient to prove theorems 1 1 . 1 1 and 1.2 in the 
situation where r\ is the identity diffeomorphism. More precisely, the F-convergence of 
ut is a consequence of the ^-convergence of (xt,yt) and the regularity of rj. Furthermore, 
from the uniform bound (depending on (xo, yo)) on the trajectories of yf) and (xt, yt) 
we deduce that g and / are uniformly bounded and Lipshitz continuous (in e, 8 < ho, 
t < min(roe,5)) over those trajectories. 

Define 

:= J 9(x,y) n(x,dy) 



9 



where \x is the family of measures introduced in Condition 1.2 Let us prove the following 
lemma. 



Lemma 7.1. 

\x nS -x nS \ <Ce Cn5 (5+{-f\+ sup \J(l)\) (7.16) 

V £ KKn 

with J(k) = Ji(k) + J 2 {k), 



k-l 



r(n+l)S 

/ g(x* n5 ,yt)ds-6g(x* nS ) (7.17) 

n=0 JnS 

and 

k-l 

:= ^2 S(g(x nS ) - g(xns, y nS )) (7.18) 

n=0 

Proof. Observe that 

r(n+l)8 r(n+l)S 

x (n+i)5 = x nS+ g{ x n&,yl)ds+ I i.9(x e s ,y e s ) - g(x e nS ,y e s ))ds (7.19) 



8 JnS 



Hence, 



with 



x (n+l)S ~ x (n+l)S = x nS ~ x nS + h + h{n) + h + h(n) + I 5 (7.20) 
r(n+l)S 

h := / (g(xl,yl)-g(x' nS ,y e s ))ds (7.21) 

JnS 

r(n+l)6 

h{n):= g(x e nS ,y e s )ds-Sg(x e nS ) (7.22) 

JnS 

h ■■= 6{g(x € nS ) - g(x nS )) (7.23) 
h(n) := 5 (g{x nS ) - g(x nS , y n s)) (7.24) 

h '■= Sg(XnS, Vns) - (X(n+1)5 ~ x ns) (7.25) 
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Now observe that 
and 



|/i| < \\Vxg\\L°°& 2 



l-^l < $\\^xg\\L°°\x e nS - X nS \. 



Using (7.7) and (7.8) we obtain that 



(7.26) 
(7.27) 

(7.28) 



Combining the previous equations, we have obtained that 

x \n+X)S ~ X (n+1)6 < x n5 ~ x nS + c{S 2 + (^f) + C5\x nS - X n5 \ + (J 2 + h)(n) (7.29) 



|/ 5 1 <c s z + 



:2 , /?\2 



and 



x \n+l)8 ~ x < 



(n+l)4 > XnS ~ x nS - C (5 2 + {-f ) - C5\x € nS - x nS \ + (I 2 + h){n) (7.30) 



Write 



n-l 



J{n) :=Y,{h + h){k) 



(7.31) 



k=0 



Summing the first n inequalities (7.29) and (7.30) 



n-l 



x € nS - x n& < c{5 2 + (^) 2 )n + C8 ^2 \x% 5 - x kS \ + J(n) 



x e n& ~ XnS > -C(8 2 + (^f)n- C8^\xU - Xks\ + J(n) 



k=0 
n-l 



Hence 



fc=0 



n-l 



\x n5 - x n5 \ < C(5 2 + ( T -f)n + C5^2\x% 5 - x k5 \ + | J(n) 

fc=0 

And we obtain by induction 

n 

K 5 - x n5 \ <C(6 2 + (-) 2 ) (n + Ctf J> - fc)(l + CS)^ 1 

6 k=l 

n 

\J{n)\+C5Y J ^ + C5) l - 2 \J{n-l + l)\ 



(7.32) 
(7.33) 

(7.34) 



(7.35) 



1=2 



Equation (7.35) concludes the proof of Lemma 7.1 



□ 



We now need to control J\{k) and J 2 (fc)- First, let us prove the following lemma. 
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Lemma 7.2. For N G N* we have 



\Ji(k)\ < (5k)C(5e c m +E 



Proof. Define y\ such that y\ = y\ for t = (n + j/N)S, j G N*, and 

#1 = 1 
(it e 

Using the regularity of / and g we obtain that 



:f«s,y e t) for (n + i/iV)5<t<(n + (j + l)/iV)<5. 



\y e t-y e t\ <C5e°Ni. 



First, observe that 



with 



and 



We have 



1 r{n+l)6 

x / ds ~ 9{x e n5 ) = Ki + K 2 

Jn5 



1 N ~ X 

E 



(n+(j+l)/N)6 



j=0 J(n+j/N)S 



{a{xns,vl) - g{x £ nS ,y e s )) ds 



iV-1 



\K-l\ < llV^Hioo— V Slip 

iV J^o (™+i/ JV )< 5 <s<("+0'+ 1 )/ w )<5 



Hence, we obtain from (7.38) that 



Moreover, we obtain from conditions 1.2 and 1.3 that 



\K 2 \ < CE( 



Ne' 



This concludes the proof of Lemma 7.2 
Lemma 7.3. We have for m G N* 



| J 2 {k)\ < C5k(m5 + E(™) + + m5 + m(^) 2 )e c ^ 
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Proof. Let m£f. Define (x 8 ,y a ) such that for j G N*, n G N*, 



= g(x s ,y s ) for jm5 < s < (j + l)m<S 
= \f(x s ,y s ) for n<5<s<n<5 + r 
y s = yn5+r for n<5 + t < s < (n + 1)5 

V{n+l)S = Vn8+r for U + 1 / jm 
^{.Xjmiyjm) = (.Xj m s,yj m s) 



Define by 



and define x° by 
Observe that 
with 



%- = lf(xjm5, Vt) for < t < (j + l)mr 

y'jniT = Vjm5 

= Xj m s for jm < n < (j + l)m. 
J 2 (A;) = K 3 + K 4 + K 5 + K 6 + K 7 

fc- 1 /-(n+l)<5 

K 3-=^2( g(x s ,y s )ds - 5g(x nS ,y n s)), 

n=0 ^ 

^=£</ 

n=0 Jn 

j> /<(n+l)T 

^:=-E( r ^n)- / g(K,ys)ds), 

T n •/ tit 



a{xn^y a s ) ds - t / g(x s ,y s )ds), 

° JnS 



k-l 



k & -.= (g(xns) -g{x a n )). 



n=0 



Using the regularity of g we obtain 



(7.46) 



(7.47) 

(7.48) 
(7.49) 

(7.50) 

(7.51) 

(7.52) 

(7.53) 



\K 6 \< SkC5m. (7.54) 
Arranging the right hand side of (7.51) into groups of m terms corresponding to the 



intervals of (7.47) we obtain, from Condition 1.2 and Condition 1.3, that 



tut 

\K 5 \ < Ck5E( — ). 



Using ( 7.48 ) and the regularity of / and g we obtain the following inequality 

\y%-m\ <Cm5e c ^. 



(7.55) 



(7.56) 
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It follows that 



\K A \ < C5kmSe cl 



Similarly, using (7.7) and (7.8), we obtain the following inequalities 



T „ 0s TXIT Qmr 



\y n s - Vns\ < C(- +m5 + m(-) ) e 



is - x n s\ < Cm(5 + (-) 2 ). 



It follows that 



\K 3 \ < CSk(^ + m5 + m(-) 2 )e 



This concludes the proof of Lemma |7.3 



(7.57) 

(7.58) 
(7.59) 

(7.60) 
□ 



Combining Lemma |7.1[ |7.2| and 7.3 we have obtained that 



\x e nS -x nS \ <Ce c5n (5+ {-f- r +5e c -k +E 



V 6 

+ (- + m5 + m(-) 2 )e- 



T 

e 

s . i 



(7.61) 



Choosing N such that e ^ ~ S 2 (observe that we need e < 5/{— Cm J)) and m such 
that ^ e c ^ ~ + we obtain for ^ + I < 1 that 

e V T e / r e — 



\xU - x n6 \ <Ce CSn (^+(-/- 6 +E(^ ^ J) 



(7.62) 



consequence of (7.2). 



This concludes the proof of inequality (7.2). The proof of (7.3) is similar and is also a 



7.2 Proof of Theorem I37T1 

Define the process t h-> (x t , y t ) by 



(x t ,y t ) := »7(itt). 



(7.63) 



It follows from the regularity of r\ that it is sufficient to prove the F-convergence of 
(xt,yt) towards 5x t <8> fJ.(X t ,dy). Now define ^ by 



Define tfj? by 



*l>h(x,y,u) : =1°C('>w) o rf l {x,y). 



(7.64) 
(7.65) 
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Proposition 7.2. The vector fields f, g and matrix fields a, Q associated with the 



system of Equations (3.2) are uniformly bounded and Lipschitz continuous. We also 
have 



(x ,yo) = r/(«o) 

XxuVt) = (xkS,yks) for kS < t < (k + 1)5 



(7.66) 



where ujk,uj' k are i.i.d. samples from the probability space (Cl,J-, P). Moreover there 
exists C > and and d-dimensional centered Gaussian vectors with identity 

covariance matrices such that for h < ho and - < tq we have 



E 



\^ 9 h( x ,V,u) - (x,y) - h(g(x,y),6) - Vh(a{x,y)£'(u),0) 



<Chi, (7.67) 



E 



\i/)t(x, y, u) - (x, y) - r(g(x, y), 0) - - (0, f(x, y)) - Vr(a(x, 0) 



[0,Q(x,y)e(u>))\ 



(7.68) 



Proof. Since (x, y) = rj(u), we obtain from (3.1) and Ito's formula 

dx = ((G + -F)Vrf o rf\x, y)) dt + (Vrf(H + —pK)) o r,~ l (x, y) dW t 
e ye 

+ \ E d ^V X {(H + ^K)(H + ^K) T ) l3 dts 
ij v v 

dy =((G + -F)Vrf o rj'^x, y)) dt + (Vr] y (H + -t=K)) o ^(x, y) dW t 
e V e 

+ ( \ E w((# + \ K )( H + -V) T L) ° t 1 dt - 



Hence we deduce from Equation (|3.2|) of Condition 3.1 that 



g(x,y) = (GVr, x + lY, d ^ X (HH T ) lJ ) o V -\x,y) 



a{x,y) = (Vt] x H) or/ 1 (x,y) 
f(x,y) = (FV V y + l^d l d jV y(KK T ) ij ) or,-\x,y) 



ij 



Q(x,y) = (Vr j y K)ori- 1 (x,y). 



(7.69) 



(7.70) 



(7.71) 

(7.72) 
(7.73) 

(7.74) 
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Remark 7.1. Observe that Equation (3.2) of Condition 3.1 requires that 

FVrf = GVr/ y , 



ij 

Y,didjrf[HH T ).. = Q 

ij 



0. 



and 



Y J did j r ] y (KH T + HK T ) =0. 



(7.75) 
(7.76) 

(7.77) 

(7.78) 

(7.79) 



ij 



Equations (7.78) and (7.79) are satisfied if KH T is skew-symmetric. One particular case 
could be, of course, KH 1 = 0, which translates into the fact that for all u the ranges of 
H{u) and K(u) are orthogonal, i.e., the noise with amplitude 1/y/e is applied to degrees 
of freedom orthogonal to those with 0(1) noise. 

We deduce the regularity of /, g, a and Q from the regularity of G, F, H, K and 
rj. Equation (7.6) is a direct consequence of the definition of ip* and ip? and Equation 
(7.66). Now observe that 

il>r(. x ,y,u) - (x,y) - T{g(x,y),0) - -(0,f(x,y)) - y/r(a(x,y)£'(u)),0) 

- J^(0, Q(x, y)?(u)) = (r,o6%-r,- r{GV V x + ^ d.drf {HH T ) l3 ,V) 



T 

e 



(7.80) 



^(0,VWH))or ? - 1 (x,y). 



Using Equations ( |7.75[ ), ( |7.76[ ), ( |7.77[ ), ( |7.78[ ) and ( |7.79[ ), the Taylor-Ito expansion of 
r/o^, the regularity of 77, and Setting £' equal to ^ defined in Equation (3.14) we obtain 
Equation (7.68). The proof of Equation (7.67) is similar. □ 



It follows from Proposition |7.2| that it is sufficient to prove Theorem 3.1 in the 
situation where rj is the identity diffeomorphism. More precisely the F-convergence of 
ut is a consequence of the F-convergence of (xt,y~t) and the regularity of r]. 

Let x 1 — y ty{x) be a function with continuous and bounded derivatives up to order 3. 
Let us prove the following lemma. 
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Lemma 7.4. We have 

E [^(^(n+l)5)] -^[¥(Xn5)] 



5E 



with 



g(xnS,yn6)V(p(x nS ) + aa T (x n s,y n s) ■ Hess ip(x n s) 



\I \<c(6l + [-) 3 > 



(7.81) 



(7.82) 



Proof. Write (x n< 5 +T , y n 5+r) := ipT( x n6,yn8,u n ). Using Equation (7.68) we obtain that 
there exists an Af(0, 1) random vector £ n independent from (x n s,y n s) and such that 



x nS+T - x nS = g{x n5 )T + y/ra(x nS , y n s)£ n + h 



with 



Hence 



(E[(/i) a ])*<C^)*. 



(7.83) 
(7.84) 



E[ip(x nS+T )] -E[ip(x n5 )] -tE g{x nS ,y n s)Vtp(x nS 
+ o-a T (x n s,y n &) ■ Hess ip(x nS ) 



(7.85) 



Similarly, using Equation (7.67) we obtain that there exists an jV(0, 1) random vector 
£' n , independent from (x n s+ T , ynS+r), and such that 



(n+l)<5 x n 



S+T = g(x nS+T , y n s +T )(5 - r) + cr(x nS+T , y„s+ T )V 6 - r£' n + I 2 (7.86) 



with 



(E[(I 2 ) 2 ]p <C(6-t)*. 



(7.87) 



Whence 



E[^(a?( n+ i) 5 )] - E[v9(x n5+r )] - (5-r)E g(x n s +T , y n s +T )Vip(x n s+ T ) 



cra T (x n s +T ,y n s +T ) : Hess ip(x n s +T ) 



(7.88) 



Using the regularity of a, we obtain that 



E 



< C(6-t) 



(7.89) 



The proof of fl7i8l) follows from ( |7!68| ), ( [7785] ), ( [7788] ), (J7789]) and the regularity of 5 and 
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Lemma 7.5. We have 



with (for 5 < Cr/e) 



n5 



Lcp(x ) 



< h 



Proof. Define Bt by -Bo = and 

B t - B nr = B n5+t - B nS for nr <t<(n + l)r. 
Define y s by y = y and 

dyt = 1 f(xo, Vt) dt + -^Q(x , y t )dBf 



Write 



g(x ) :=J g(x ,y) n{xQ,dy). 



Using Lemma 7.4 we obtain 

E[</?(x n(5 )] - ip(x ) 



n5 



Lif(x ) + Ji + J 2 + J 3 + Ja, 



with 



L(f(x ) := g{x )Vtp(x ) + cra T (xo) : Ressip(x ), 

Ji =- Y] E g(x k s, yks)^^p(xks) + o-o- T (x kS , y k s) ■ Hess <p(x k s) 
n L 
k=o 

^ n— 1 

L(xo,2/fc,5)V^(xo) + aa T (x Q ,y k s) ■ Hess</?(x ) 

fc=o 

^ n— 1 

^ 2 = n ^ ( E [9 , (^o,yfc5)Vv9(x ) + aa T {xo,y kS ) : Hess^(xo) 

5(x ,y s )V^(j;o) + o-0- T Oro,2/s) : Hess<^(x ) 



fc=o 

I r(k+l)r 



E 



r ifcr 



1 f nT 

Ji = — E 



(is 



nr Jo 



g(x , y s )Vip(x ) + aa T (x , y s ) ■ Hess <^(x ) ds - L(p(x ) 



/Tx a 1 



\M < c 55 + 
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Using the regularity of a, #,99, (7.6) and (7.7) we obtain 



/ 1 T 

Ji\ < Ci(n5)2 +n5 + n{-) 



Using Property 3 of Condition 3.1 and Property 3 of Condition 3.2 we obtain 

TIT 

\M<CE{-). 



Using (7.67) and (7.68), we obtain the following inequality 

< c(J^+ {n§)% +nS + n(-) 



E 



UnS Dm 



TIT n nr 



which leads to 
Hence, we have obtained 

E[^(x„5)] - ip(x 



n5 



Ltp(x Q ) 



< Js, 



with 



I h\ < C(\l~ + {nS)\ + nd + n(~) ^e cn r + E {™) + C(~) ^. 



Choosing n such that 



— e e ~ 



we obtain (7.91) for 5 < Cr/e. 



(7.101) 



(7.102) 



(7.103) 



(7.104) 



(7.105) 



(7.106) 



□ 



We now combine Lemma |7.5| with Theorem 1 of Chapter 2 of [104] which states 
that the uniform convergence (in xq, yo) of ^ "^J — ^ — to Lip(xo) as e ! 0, r < 5, 

3 

f | 0, f I and (§) 5 J 1 implies the convergence in distribution of x n s to the Markov 
process generated by L. 

The F-convergence of (xt,yt) can be deduced from the convergence in distribution 
of xt an d Equation (3.4) of Condition 3.1 The proof follows the same lines as above, 
which will not be repeated here. 
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